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 <optional>
9 :
10 : #include "DataStructures/CachedTempBuffer.hpp"
11 : #include "DataStructures/DataBox/Prefixes.hpp"
12 : #include "DataStructures/TaggedTuple.hpp"
13 : #include "DataStructures/TempBuffer.hpp"
14 : #include "DataStructures/Tensor/EagerMath/Magnitude.hpp"
15 : #include "DataStructures/Tensor/Tensor.hpp"
16 : #include "Elliptic/Systems/Xcts/Tags.hpp"
17 : #include "NumericalAlgorithms/LinearOperators/PartialDerivatives.hpp"
18 : #include "Options/Auto.hpp"
19 : #include "Options/Context.hpp"
20 : #include "Options/ParseError.hpp"
21 : #include "Options/String.hpp"
22 : #include "PointwiseFunctions/AnalyticData/Xcts/CommonVariables.hpp"
23 : #include "PointwiseFunctions/AnalyticSolutions/Xcts/Flatness.hpp"
24 : #include "PointwiseFunctions/GeneralRelativity/Tags.hpp"
25 : #include "PointwiseFunctions/GeneralRelativity/Tags/Conformal.hpp"
26 : #include "PointwiseFunctions/InitialDataUtilities/Background.hpp"
27 : #include "PointwiseFunctions/InitialDataUtilities/InitialGuess.hpp"
28 : #include "Utilities/CallWithDynamicType.hpp"
29 : #include "Utilities/Requires.hpp"
30 : #include "Utilities/Serialization/CharmPupable.hpp"
31 : #include "Utilities/Serialization/PupStlCpp17.hpp"
32 : #include "Utilities/TMPL.hpp"
33 :
34 : /// \cond
35 : namespace PUP {
36 : class er;
37 : } // namespace PUP
38 : /// \endcond
39 :
40 : namespace Xcts::AnalyticData {
41 :
42 : namespace detail {
43 :
44 : template <typename DataType>
45 : using BinaryVariablesCache = cached_temp_buffer_from_typelist<tmpl::append<
46 : common_tags<DataType>,
47 : tmpl::list<gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
48 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
49 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, 3>, 0>,
50 : // For initial guesses
51 : Tags::ConformalFactorMinusOne<DataType>,
52 : Tags::LapseTimesConformalFactorMinusOne<DataType>,
53 : Tags::ShiftExcess<DataType, 3, Frame::Inertial>>,
54 : hydro_tags<DataType>>>;
55 :
56 : template <typename DataType>
57 : struct BinaryVariables
58 : : CommonVariables<DataType, BinaryVariablesCache<DataType>> {
59 : static constexpr size_t Dim = 3;
60 : using Cache = BinaryVariablesCache<DataType>;
61 : using Base = CommonVariables<DataType, BinaryVariablesCache<DataType>>;
62 : using Base::operator();
63 :
64 : using superposed_tags = tmpl::append<
65 : tmpl::list<
66 : Tags::ConformalMetric<DataType, Dim, Frame::Inertial>,
67 : ::Tags::deriv<Tags::ConformalMetric<DataType, Dim, Frame::Inertial>,
68 : tmpl::size_t<Dim>, Frame::Inertial>,
69 : gr::Tags::TraceExtrinsicCurvature<DataType>,
70 : ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataType>>,
71 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0>,
72 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0>,
73 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, Dim>, 0>,
74 : Tags::ConformalFactorMinusOne<DataType>,
75 : Tags::LapseTimesConformalFactorMinusOne<DataType>,
76 : Tags::ShiftExcess<DataType, Dim, Frame::Inertial>>,
77 : hydro_tags<DataType>>;
78 :
79 : BinaryVariables(
80 : std::optional<std::reference_wrapper<const Mesh<Dim>>> local_mesh,
81 : std::optional<std::reference_wrapper<const InverseJacobian<
82 : DataType, Dim, Frame::ElementLogical, Frame::Inertial>>>
83 : local_inv_jacobian,
84 : const tnsr::I<DataVector, Dim>& local_x,
85 : const double local_angular_velocity, const double local_expansion,
86 : const std::array<double, 3> local_linear_velocity,
87 : std::optional<std::array<double, 2>> local_falloff_widths,
88 : std::array<tnsr::I<DataVector, Dim>, 2> local_x_isolated,
89 : std::array<DataVector, 2> local_windows,
90 : tuples::tagged_tuple_from_typelist<superposed_tags> local_flat_vars,
91 : std::array<tuples::tagged_tuple_from_typelist<superposed_tags>, 2>
92 : local_isolated_vars)
93 : : Base(std::move(local_mesh), std::move(local_inv_jacobian)),
94 : x(local_x),
95 : angular_velocity(local_angular_velocity),
96 : expansion(local_expansion),
97 : linear_velocity(local_linear_velocity),
98 : falloff_widths(std::move(local_falloff_widths)),
99 : x_isolated(std::move(local_x_isolated)),
100 : windows(std::move(local_windows)),
101 : flat_vars(std::move(local_flat_vars)),
102 : isolated_vars(std::move(local_isolated_vars)) {}
103 :
104 : const tnsr::I<DataVector, Dim>& x;
105 : const double angular_velocity;
106 : const double expansion;
107 : const std::array<double, 3> linear_velocity;
108 : const std::optional<std::array<double, 2>> falloff_widths;
109 : const std::array<tnsr::I<DataVector, Dim>, 2> x_isolated;
110 : const std::array<DataVector, 2> windows;
111 : const tuples::tagged_tuple_from_typelist<superposed_tags> flat_vars;
112 : const std::array<tuples::tagged_tuple_from_typelist<superposed_tags>, 2>
113 : isolated_vars;
114 :
115 : template <bool ApplyWindow = true, typename Tag,
116 : Requires<tmpl::list_contains_v<superposed_tags, Tag>> = nullptr>
117 : void superposition(gsl::not_null<typename Tag::type*> superposed_var,
118 : gsl::not_null<Cache*> /*cache*/, Tag /*meta*/) const {
119 : for (size_t i = 0; i < superposed_var->size(); ++i) {
120 : if constexpr (ApplyWindow) {
121 : (*superposed_var)[i] =
122 : get<Tag>(flat_vars)[i] +
123 : windows[0] *
124 : (get<Tag>(isolated_vars[0])[i] - get<Tag>(flat_vars)[i]) +
125 : windows[1] *
126 : (get<Tag>(isolated_vars[1])[i] - get<Tag>(flat_vars)[i]);
127 : } else {
128 : (*superposed_var)[i] = get<Tag>(isolated_vars[0])[i] +
129 : get<Tag>(isolated_vars[1])[i] -
130 : get<Tag>(flat_vars)[i];
131 : }
132 : }
133 : }
134 :
135 : void operator()(
136 : const gsl::not_null<tnsr::ii<DataType, Dim>*> conformal_metric,
137 : const gsl::not_null<Cache*> cache,
138 : Tags::ConformalMetric<DataType, Dim, Frame::Inertial> meta)
139 : const override {
140 : superposition(conformal_metric, cache, meta);
141 : }
142 : void operator()(
143 : const gsl::not_null<tnsr::ijj<DataType, Dim>*> deriv_conformal_metric,
144 : const gsl::not_null<Cache*> cache,
145 : ::Tags::deriv<Tags::ConformalMetric<DataType, Dim, Frame::Inertial>,
146 : tmpl::size_t<Dim>, Frame::Inertial>
147 : meta) const override {
148 : superposition(deriv_conformal_metric, cache, meta);
149 : add_deriv_of_window_function(deriv_conformal_metric);
150 : }
151 : void operator()(
152 : const gsl::not_null<Scalar<DataType>*> extrinsic_curvature_trace,
153 : const gsl::not_null<Cache*> cache,
154 : gr::Tags::TraceExtrinsicCurvature<DataType> meta) const override {
155 : superposition(extrinsic_curvature_trace, cache, meta);
156 : }
157 : void operator()(
158 : const gsl::not_null<Scalar<DataType>*> dt_extrinsic_curvature_trace,
159 : const gsl::not_null<Cache*> cache,
160 : ::Tags::dt<gr::Tags::TraceExtrinsicCurvature<DataType>> meta)
161 : const override {
162 : superposition(dt_extrinsic_curvature_trace, cache, meta);
163 : }
164 : void operator()(
165 : gsl::not_null<tnsr::I<DataType, Dim>*> shift_background,
166 : gsl::not_null<Cache*> cache,
167 : Tags::ShiftBackground<DataType, Dim, Frame::Inertial> /*meta*/)
168 : const override;
169 : void operator()(
170 : gsl::not_null<tnsr::iJ<DataType, Dim>*> deriv_shift_background,
171 : gsl::not_null<Cache*> cache,
172 : ::Tags::deriv<Tags::ShiftBackground<DataType, Dim, Frame::Inertial>,
173 : tmpl::size_t<Dim>, Frame::Inertial> /*meta*/)
174 : const override;
175 : void operator()(
176 : const gsl::not_null<Scalar<DataType>*> conformal_energy_density,
177 : const gsl::not_null<Cache*> cache,
178 : gr::Tags::Conformal<gr::Tags::EnergyDensity<DataType>, 0> meta) const {
179 : superposition<false>(conformal_energy_density, cache, meta);
180 : }
181 : void operator()(
182 : const gsl::not_null<Scalar<DataType>*> conformal_stress_trace,
183 : const gsl::not_null<Cache*> cache,
184 : gr::Tags::Conformal<gr::Tags::StressTrace<DataType>, 0> meta) const {
185 : superposition<false>(conformal_stress_trace, cache, meta);
186 : }
187 : void operator()(
188 : const gsl::not_null<tnsr::I<DataType, Dim>*> conformal_momentum_density,
189 : const gsl::not_null<Cache*> cache,
190 : gr::Tags::Conformal<gr::Tags::MomentumDensity<DataType, Dim>, 0> meta)
191 : const {
192 : superposition<false>(conformal_momentum_density, cache, meta);
193 : }
194 : void operator()(
195 : const gsl::not_null<Scalar<DataType>*> conformal_factor_minus_one,
196 : const gsl::not_null<Cache*> cache,
197 : Tags::ConformalFactorMinusOne<DataType> meta) const {
198 : superposition(conformal_factor_minus_one, cache, meta);
199 : }
200 : void operator()(
201 : const gsl::not_null<Scalar<DataType>*>
202 : lapse_times_conformal_factor_minus_one,
203 : const gsl::not_null<Cache*> cache,
204 : Tags::LapseTimesConformalFactorMinusOne<DataType> meta) const {
205 : superposition(lapse_times_conformal_factor_minus_one, cache, meta);
206 : }
207 : void operator()(
208 : const gsl::not_null<tnsr::I<DataType, Dim>*> shift_excess,
209 : const gsl::not_null<Cache*> cache,
210 : Tags::ShiftExcess<DataType, Dim, Frame::Inertial> meta) const {
211 : superposition(shift_excess, cache, meta);
212 : }
213 : void operator()(const gsl::not_null<Scalar<DataType>*> rest_mass_density,
214 : const gsl::not_null<Cache*> cache,
215 : hydro::Tags::RestMassDensity<DataType> meta) const {
216 : superposition<false>(rest_mass_density, cache, meta);
217 : }
218 : void operator()(const gsl::not_null<Scalar<DataType>*> specific_enthalpy,
219 : const gsl::not_null<Cache*> cache,
220 : hydro::Tags::SpecificEnthalpy<DataType> meta) const {
221 : superposition<false>(specific_enthalpy, cache, meta);
222 : }
223 : void operator()(const gsl::not_null<Scalar<DataType>*> pressure,
224 : const gsl::not_null<Cache*> cache,
225 : hydro::Tags::Pressure<DataType> meta) const {
226 : superposition<false>(pressure, cache, meta);
227 : }
228 : void operator()(const gsl::not_null<tnsr::I<DataType, 3>*> spatial_velocity,
229 : const gsl::not_null<Cache*> cache,
230 : hydro::Tags::SpatialVelocity<DataType, 3> meta) const {
231 : superposition<false>(spatial_velocity, cache, meta);
232 : }
233 : void operator()(const gsl::not_null<Scalar<DataType>*> lorentz_factor,
234 : const gsl::not_null<Cache*> cache,
235 : hydro::Tags::LorentzFactor<DataType> meta) const {
236 : superposition<false>(lorentz_factor, cache, meta);
237 : }
238 : void operator()(const gsl::not_null<tnsr::I<DataType, 3>*> magnetic_field,
239 : const gsl::not_null<Cache*> cache,
240 : hydro::Tags::MagneticField<DataType, 3> meta) const {
241 : superposition<false>(magnetic_field, cache, meta);
242 : }
243 :
244 : private:
245 : void add_deriv_of_window_function(
246 : gsl::not_null<tnsr::ijj<DataType, Dim>*> deriv_conformal_metric) const;
247 : };
248 : } // namespace detail
249 :
250 : /*!
251 : * \brief Binary compact-object data in general relativity, constructed from
252 : * superpositions of two isolated objects.
253 : *
254 : * This class implements background data for the XCTS equations describing two
255 : * objects in a quasi-equilibrium orbit, i.e. with \f$\bar{u}=0\f$ and
256 : * \f$\partial_t K=0\f$. Both objects can be chosen from the list of
257 : * `IsolatedObjectRegistrars`, e.g. they can be black-hole or neutron-star
258 : * solutions in different coordinates. Most quantities are constructed by
259 : * superposing the two isolated solutions (see e.g. Eq. (8-9) in
260 : * \cite Varma2018sqd or Eq. (45-46) in \cite Lovelace2008tw):
261 : *
262 : * \f{align}
263 : * \bar{\gamma}_{ij} &= f_{ij} + \sum_{\alpha=1}^2
264 : * e^{-r_\alpha^2 / w_\alpha^2}\left(\gamma^\alpha_{ij} - f_{ij}\right) \\
265 : * K &= \sum_{\alpha=1}^2 e^{-r_\alpha^2 / w_\alpha^2}K^\alpha
266 : * \f}
267 : *
268 : * where \f$\gamma^\alpha_{ij}\f$ and \f$K^\alpha\f$ denote the spatial metric
269 : * and extrinsic-curvature trace of the two individual solutions, \f$r_\alpha\f$
270 : * is the Euclidean coordinate-distance from the center of each object and
271 : * \f$w_\alpha\f$ are parameters describing the falloff widths of Gaussian
272 : * window functions. The window functions
273 : * facilitate that the influence of either of the two objects
274 : * at the position of the other is strongly damped, and they also avoid
275 : * logarithmic scaling of the solution at large distances where we would
276 : * typically employ an inverse-radial coordinate map and asymptotically-flat
277 : * boundary conditions. The falloff-widths are chosen in terms of the Newtonian
278 : * Lagrange points of the two objects in \cite Varma2018sqd and
279 : * \cite Lovelace2008tw, and they are input parameters in this implementation.
280 : * The falloff can be disabled by passing `std::nullopt` to the constructor, or
281 : * `None` in the input file.
282 : *
283 : * \par Matter sources
284 : * Matter sources are superposed without the window functions. The analytic
285 : * matter sources are of
286 : * limited use anyway, because in a binary setting they don't take the
287 : * gravitational influence of the other body into account. Therefore, the matter
288 : * sources should typically be solved-for alongside the gravity sector to impose
289 : * conditions such as hydrostatic equilibrium. For scenarios where we just want
290 : * to superpose the isolated matter solutions and compute the resulting gravity,
291 : * the matter sources are simply added.
292 : *
293 : * \par Orbital motion
294 : * The remaining quantities that this class implements relate to the orbital
295 : * motion of the two objects. To obtain initial data in "co-rotating"
296 : * coordinates where the two objects are initially at rest we prescribe the
297 : * background shift
298 : *
299 : * \f{equation} \beta^i_\mathrm{background} = (-\Omega y, \Omega x, 0) +
300 : * \dot{a}_0 x^i + v^i_0 \f}
301 : *
302 : * where \f$\Omega\f$ is the angular-velocity parameter and \f$\dot{a}_0\f$
303 : * is an expansion parameter. Both control the eccentricity of the orbit.
304 : * The parameter \f$v^i_0\f$ is a constant velocity that can be used to
305 : * control the linear momentum of the system (see Eq. (28) in
306 : * \cite Ossokine2015yla).
307 : */
308 : template <typename IsolatedObjectBase, typename IsolatedObjectClasses>
309 1 : class Binary : public elliptic::analytic_data::Background,
310 : public elliptic::analytic_data::InitialGuess {
311 : public:
312 : template <typename DataType>
313 0 : using tags = typename detail::BinaryVariablesCache<DataType>::tags_list;
314 :
315 0 : struct XCoords {
316 0 : static constexpr Options::String help =
317 : "The coordinates on the x-axis where the two objects are placed";
318 0 : using type = std::array<double, 2>;
319 : };
320 0 : struct CenterOfMassOffset {
321 0 : static constexpr Options::String help = {
322 : "Offset in the y and z axes applied to both objects in order to "
323 : "control the center of mass."};
324 0 : using type = std::array<double, 2>;
325 : };
326 0 : struct ObjectLeft {
327 0 : static constexpr Options::String help =
328 : "The object placed on the negative x-axis";
329 0 : using type = std::unique_ptr<IsolatedObjectBase>;
330 : };
331 0 : struct ObjectRight {
332 0 : static constexpr Options::String help =
333 : "The object placed on the positive x-axis";
334 0 : using type = std::unique_ptr<IsolatedObjectBase>;
335 : };
336 0 : struct AngularVelocity {
337 0 : static constexpr Options::String help =
338 : "Orbital angular velocity 'Omega0' about the z-axis. Added to the "
339 : "background shift as a term 'Omega0 x r'.";
340 0 : using type = double;
341 : };
342 0 : struct Expansion {
343 0 : static constexpr Options::String help =
344 : "The expansion parameter 'adot0', which is a radial velocity over "
345 : "radius. Added to the background shift as a term 'adot0 r^i'";
346 0 : using type = double;
347 : };
348 0 : struct LinearVelocity {
349 0 : static constexpr Options::String help =
350 : "Constant velocity 'v0' added to the background shift to control the "
351 : "linear momentum of the system.";
352 0 : using type = std::array<double, 3>;
353 : };
354 0 : struct FalloffWidths {
355 0 : static constexpr Options::String help =
356 : "The widths for the window functions around the two objects, or 'None' "
357 : "to disable the Gaussian falloff.";
358 0 : using type = Options::Auto<std::array<double, 2>, Options::AutoLabel::None>;
359 : };
360 0 : using options =
361 : tmpl::list<XCoords, CenterOfMassOffset, ObjectLeft, ObjectRight,
362 : AngularVelocity, Expansion, LinearVelocity, FalloffWidths>;
363 0 : static constexpr Options::String help =
364 : "Binary compact-object data in general relativity, constructed from "
365 : "superpositions of two isolated objects.";
366 :
367 0 : Binary() = default;
368 0 : Binary(const Binary&) = delete;
369 0 : Binary& operator=(const Binary&) = delete;
370 0 : Binary(Binary&&) = default;
371 0 : Binary& operator=(Binary&&) = default;
372 0 : ~Binary() = default;
373 :
374 0 : Binary(const std::array<double, 2> xcoords,
375 : const std::array<double, 2> center_of_mass_offset,
376 : std::unique_ptr<IsolatedObjectBase> object_left,
377 : std::unique_ptr<IsolatedObjectBase> object_right,
378 : const double angular_velocity, const double expansion,
379 : const std::array<double, 3> linear_velocity,
380 : const std::optional<std::array<double, 2>> falloff_widths,
381 : const Options::Context& context = {})
382 : : xcoords_(xcoords),
383 : y_offset_(center_of_mass_offset[0]),
384 : z_offset_(center_of_mass_offset[1]),
385 : superposed_objects_({std::move(object_left), std::move(object_right)}),
386 : angular_velocity_(angular_velocity),
387 : expansion_(expansion),
388 : linear_velocity_(linear_velocity),
389 : falloff_widths_(falloff_widths) {
390 : if (xcoords_[0] >= xcoords_[1]) {
391 : PARSE_ERROR(context, "Specify 'XCoords' ascending from left to right.");
392 : }
393 : }
394 :
395 0 : explicit Binary(CkMigrateMessage* m)
396 : : elliptic::analytic_data::Background(m),
397 : elliptic::analytic_data::InitialGuess(m) {}
398 : using PUP::able::register_constructor;
399 0 : WRAPPED_PUPable_decl_template(Binary);
400 :
401 : template <typename DataType, typename... RequestedTags,
402 : Requires<tmpl2::flat_all_v<tmpl::list_contains_v<
403 : tags<DataType>, RequestedTags>...>> = nullptr>
404 0 : tuples::TaggedTuple<RequestedTags...> variables(
405 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
406 : tmpl::list<RequestedTags...> /*meta*/) const {
407 : return variables_impl<DataType>(x, std::nullopt, std::nullopt,
408 : tmpl::list<RequestedTags...>{});
409 : }
410 : template <typename... RequestedTags,
411 : Requires<tmpl2::flat_all_v<tmpl::list_contains_v<
412 : tags<DataVector>, RequestedTags>...>> = nullptr>
413 0 : tuples::TaggedTuple<RequestedTags...> variables(
414 : const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
415 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
416 : Frame::Inertial>& inv_jacobian,
417 : tmpl::list<RequestedTags...> /*meta*/) const {
418 : return variables_impl<DataVector>(x, mesh, inv_jacobian,
419 : tmpl::list<RequestedTags...>{});
420 : }
421 :
422 : // NOLINTNEXTLINE
423 0 : void pup(PUP::er& p) override {
424 : elliptic::analytic_data::Background::pup(p);
425 : elliptic::analytic_data::InitialGuess::pup(p);
426 : p | xcoords_;
427 : p | y_offset_;
428 : p | z_offset_;
429 : p | superposed_objects_;
430 : p | angular_velocity_;
431 : p | expansion_;
432 : p | linear_velocity_;
433 : p | falloff_widths_;
434 : }
435 :
436 : /// Coordinates of the objects, ascending left to right
437 1 : const std::array<double, 2>& x_coords() const { return xcoords_; }
438 : /// Offset in y and z coordinates of the objects
439 1 : double y_offset() const { return y_offset_; }
440 0 : double z_offset() const { return z_offset_; }
441 : /// The two objects. First entry is the left object, second entry is the right
442 : /// object.
443 1 : const std::array<std::unique_ptr<IsolatedObjectBase>, 2>& superposed_objects()
444 : const {
445 : return superposed_objects_;
446 : }
447 0 : double angular_velocity() const { return angular_velocity_; }
448 0 : double expansion() const { return expansion_; }
449 0 : const std::array<double, 3>& linear_velocity() const {
450 : return linear_velocity_;
451 : }
452 0 : const std::optional<std::array<double, 2>>& falloff_widths() const {
453 : return falloff_widths_;
454 : }
455 :
456 : private:
457 0 : std::array<double, 2> xcoords_{};
458 0 : double y_offset_{};
459 0 : double z_offset_{};
460 0 : std::array<std::unique_ptr<IsolatedObjectBase>, 2> superposed_objects_{};
461 0 : Xcts::Solutions::Flatness flatness_{};
462 0 : double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
463 0 : double expansion_ = std::numeric_limits<double>::signaling_NaN();
464 0 : std::array<double, 3> linear_velocity_{};
465 0 : std::optional<std::array<double, 2>> falloff_widths_{};
466 :
467 : template <typename DataType, typename... RequestedTags,
468 : Requires<tmpl2::flat_all_v<tmpl::list_contains_v<
469 : tags<DataType>, RequestedTags>...>> = nullptr>
470 0 : tuples::TaggedTuple<RequestedTags...> variables_impl(
471 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
472 : std::optional<std::reference_wrapper<const Mesh<3>>> mesh,
473 : std::optional<std::reference_wrapper<const InverseJacobian<
474 : DataType, 3, Frame::ElementLogical, Frame::Inertial>>>
475 : inv_jacobian,
476 : tmpl::list<RequestedTags...> /*meta*/) const {
477 : std::array<tnsr::I<DataVector, 3>, 2> x_isolated{{x, x}};
478 : const std::array<std::array<double, 3>, 2> coords_isolated{
479 : {{{xcoords_[0], y_offset_, z_offset_}},
480 : {{xcoords_[1], y_offset_, z_offset_}}}};
481 : std::array<DataVector, 2> euclidean_distance{};
482 : std::array<DataVector, 2> windows{};
483 : // Possible optimization: Only retrieve those superposed tags from the
484 : // isolated solutions that are actually needed. This needs some dependency
485 : // logic, because some of the non-superposed tags depend on superposed tags.
486 : using VarsComputer = detail::BinaryVariables<DataType>;
487 : using requested_superposed_tags = typename VarsComputer::superposed_tags;
488 : std::array<tuples::tagged_tuple_from_typelist<requested_superposed_tags>, 2>
489 : isolated_vars;
490 : for (size_t i = 0; i < 2; ++i) {
491 : for (size_t dim = 0; dim < 3; dim++) {
492 : gsl::at(x_isolated, i).get(dim) -= gsl::at(coords_isolated, i)[dim];
493 : }
494 : gsl::at(euclidean_distance, i) = get(magnitude(gsl::at(x_isolated, i)));
495 : if (falloff_widths_.has_value()) {
496 : gsl::at(windows, i) = exp(-square(gsl::at(euclidean_distance, i)) /
497 : square(gsl::at(*falloff_widths_, i)));
498 : } else {
499 : gsl::at(windows, i) = make_with_value<DataVector>(x, 1.);
500 : }
501 : gsl::at(isolated_vars, i) = get_isolated_vars<requested_superposed_tags>(
502 : *gsl::at(superposed_objects_, i), gsl::at(x_isolated, i));
503 : }
504 : auto flat_vars = flatness_.variables(x, requested_superposed_tags{});
505 : typename VarsComputer::Cache cache{get_size(*x.begin())};
506 : const VarsComputer computer{std::move(mesh),
507 : std::move(inv_jacobian),
508 : x,
509 : angular_velocity_,
510 : expansion_,
511 : linear_velocity_,
512 : falloff_widths_,
513 : std::move(x_isolated),
514 : std::move(windows),
515 : std::move(flat_vars),
516 : std::move(isolated_vars)};
517 : return {cache.get_var(computer, RequestedTags{})...};
518 : }
519 :
520 : template <typename TagsList, typename... Args>
521 0 : tuples::tagged_tuple_from_typelist<TagsList> get_isolated_vars(
522 : const IsolatedObjectBase& isolated_object, const Args&... args) const {
523 : return call_with_dynamic_type<tuples::tagged_tuple_from_typelist<TagsList>,
524 : IsolatedObjectClasses>(
525 : &isolated_object, [&args...](const auto* const derived) {
526 : return derived->variables(args..., TagsList{});
527 : });
528 : }
529 : };
530 :
531 : /// \cond
532 : template <typename IsolatedObjectBase, typename IsolatedObjectClasses>
533 : PUP::able::PUP_ID Binary<IsolatedObjectBase, IsolatedObjectClasses>::my_PUP_ID =
534 : 0; // NOLINT
535 : /// \endcond
536 :
537 : } // namespace Xcts::AnalyticData
|
