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