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<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 0 : struct XCoords {
313 0 : static constexpr Options::String help =
314 : "The coordinates on the x-axis where the two objects are placed";
315 0 : using type = std::array<double, 2>;
316 : };
317 0 : struct CenterOfMassOffset {
318 0 : static constexpr Options::String help = {
319 : "Offset in the y and z axes applied to both objects in order to "
320 : "control the center of mass."};
321 0 : using type = std::array<double, 2>;
322 : };
323 0 : struct ObjectLeft {
324 0 : static constexpr Options::String help =
325 : "The object placed on the negative x-axis";
326 0 : using type = std::unique_ptr<IsolatedObjectBase>;
327 : };
328 0 : struct ObjectRight {
329 0 : static constexpr Options::String help =
330 : "The object placed on the positive x-axis";
331 0 : using type = std::unique_ptr<IsolatedObjectBase>;
332 : };
333 0 : struct AngularVelocity {
334 0 : static constexpr Options::String help =
335 : "Orbital angular velocity 'Omega0' about the z-axis. Added to the "
336 : "background shift as a term 'Omega0 x r'.";
337 0 : using type = double;
338 : };
339 0 : struct Expansion {
340 0 : static constexpr Options::String help =
341 : "The expansion parameter 'adot0', which is a radial velocity over "
342 : "radius. Added to the background shift as a term 'adot0 r^i'";
343 0 : using type = double;
344 : };
345 0 : struct LinearVelocity {
346 0 : static constexpr Options::String help =
347 : "Constant velocity 'v0' added to the background shift to control the "
348 : "linear momentum of the system.";
349 0 : using type = std::array<double, 3>;
350 : };
351 0 : struct FalloffWidths {
352 0 : static constexpr Options::String help =
353 : "The widths for the window functions around the two objects, or 'None' "
354 : "to disable the Gaussian falloff.";
355 0 : using type = Options::Auto<std::array<double, 2>, Options::AutoLabel::None>;
356 : };
357 0 : using options =
358 : tmpl::list<XCoords, CenterOfMassOffset, ObjectLeft, ObjectRight,
359 : AngularVelocity, Expansion, LinearVelocity, FalloffWidths>;
360 0 : static constexpr Options::String help =
361 : "Binary compact-object data in general relativity, constructed from "
362 : "superpositions of two isolated objects.";
363 :
364 0 : Binary() = default;
365 0 : Binary(const Binary&) = delete;
366 0 : Binary& operator=(const Binary&) = delete;
367 0 : Binary(Binary&&) = default;
368 0 : Binary& operator=(Binary&&) = default;
369 0 : ~Binary() = default;
370 :
371 0 : Binary(const std::array<double, 2> xcoords,
372 : const std::array<double, 2> center_of_mass_offset,
373 : std::unique_ptr<IsolatedObjectBase> object_left,
374 : std::unique_ptr<IsolatedObjectBase> object_right,
375 : const double angular_velocity, const double expansion,
376 : const std::array<double, 3> linear_velocity,
377 : const std::optional<std::array<double, 2>> falloff_widths,
378 : const Options::Context& context = {})
379 : : xcoords_(xcoords),
380 : y_offset_(center_of_mass_offset[0]),
381 : z_offset_(center_of_mass_offset[1]),
382 : superposed_objects_({std::move(object_left), std::move(object_right)}),
383 : angular_velocity_(angular_velocity),
384 : expansion_(expansion),
385 : linear_velocity_(linear_velocity),
386 : falloff_widths_(falloff_widths) {
387 : if (xcoords_[0] >= xcoords_[1]) {
388 : PARSE_ERROR(context, "Specify 'XCoords' ascending from left to right.");
389 : }
390 : }
391 :
392 0 : explicit Binary(CkMigrateMessage* m)
393 : : elliptic::analytic_data::Background(m),
394 : elliptic::analytic_data::InitialGuess(m) {}
395 : using PUP::able::register_constructor;
396 0 : WRAPPED_PUPable_decl_template(Binary);
397 :
398 : template <typename DataType, typename... RequestedTags>
399 0 : tuples::TaggedTuple<RequestedTags...> variables(
400 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
401 : tmpl::list<RequestedTags...> /*meta*/) const {
402 : return variables_impl<DataType>(x, std::nullopt, std::nullopt,
403 : tmpl::list<RequestedTags...>{});
404 : }
405 : template <typename... RequestedTags>
406 0 : tuples::TaggedTuple<RequestedTags...> variables(
407 : const tnsr::I<DataVector, 3, Frame::Inertial>& x, const Mesh<3>& mesh,
408 : const InverseJacobian<DataVector, 3, Frame::ElementLogical,
409 : Frame::Inertial>& inv_jacobian,
410 : tmpl::list<RequestedTags...> /*meta*/) const {
411 : return variables_impl<DataVector>(x, mesh, inv_jacobian,
412 : tmpl::list<RequestedTags...>{});
413 : }
414 :
415 : // NOLINTNEXTLINE
416 0 : void pup(PUP::er& p) override {
417 : elliptic::analytic_data::Background::pup(p);
418 : elliptic::analytic_data::InitialGuess::pup(p);
419 : p | xcoords_;
420 : p | y_offset_;
421 : p | z_offset_;
422 : p | superposed_objects_;
423 : p | angular_velocity_;
424 : p | expansion_;
425 : p | linear_velocity_;
426 : p | falloff_widths_;
427 : }
428 :
429 : /// Coordinates of the objects, ascending left to right
430 1 : const std::array<double, 2>& x_coords() const { return xcoords_; }
431 : /// Offset in y and z coordinates of the objects
432 1 : double y_offset() const { return y_offset_; }
433 0 : double z_offset() const { return z_offset_; }
434 : /// The two objects. First entry is the left object, second entry is the right
435 : /// object.
436 1 : const std::array<std::unique_ptr<IsolatedObjectBase>, 2>& superposed_objects()
437 : const {
438 : return superposed_objects_;
439 : }
440 0 : double angular_velocity() const { return angular_velocity_; }
441 0 : double expansion() const { return expansion_; }
442 0 : const std::array<double, 3>& linear_velocity() const {
443 : return linear_velocity_;
444 : }
445 0 : const std::optional<std::array<double, 2>>& falloff_widths() const {
446 : return falloff_widths_;
447 : }
448 :
449 : private:
450 0 : std::array<double, 2> xcoords_{};
451 0 : double y_offset_{};
452 0 : double z_offset_{};
453 0 : std::array<std::unique_ptr<IsolatedObjectBase>, 2> superposed_objects_{};
454 0 : Xcts::Solutions::Flatness flatness_{};
455 0 : double angular_velocity_ = std::numeric_limits<double>::signaling_NaN();
456 0 : double expansion_ = std::numeric_limits<double>::signaling_NaN();
457 0 : std::array<double, 3> linear_velocity_{};
458 0 : std::optional<std::array<double, 2>> falloff_widths_{};
459 :
460 : template <typename DataType, typename... RequestedTags>
461 0 : tuples::TaggedTuple<RequestedTags...> variables_impl(
462 : const tnsr::I<DataType, 3, Frame::Inertial>& x,
463 : std::optional<std::reference_wrapper<const Mesh<3>>> mesh,
464 : std::optional<std::reference_wrapper<const InverseJacobian<
465 : DataType, 3, Frame::ElementLogical, Frame::Inertial>>>
466 : inv_jacobian,
467 : tmpl::list<RequestedTags...> /*meta*/) const {
468 : std::array<tnsr::I<DataVector, 3>, 2> x_isolated{{x, x}};
469 : const std::array<std::array<double, 3>, 2> coords_isolated{
470 : {{{xcoords_[0], y_offset_, z_offset_}},
471 : {{xcoords_[1], y_offset_, z_offset_}}}};
472 : std::array<DataVector, 2> euclidean_distance{};
473 : std::array<DataVector, 2> windows{};
474 : // Possible optimization: Only retrieve those superposed tags from the
475 : // isolated solutions that are actually needed. This needs some dependency
476 : // logic, because some of the non-superposed tags depend on superposed tags.
477 : using VarsComputer = detail::BinaryVariables<DataType>;
478 : using requested_superposed_tags = typename VarsComputer::superposed_tags;
479 : std::array<tuples::tagged_tuple_from_typelist<requested_superposed_tags>, 2>
480 : isolated_vars;
481 : for (size_t i = 0; i < 2; ++i) {
482 : for (size_t dim = 0; dim < 3; dim++) {
483 : gsl::at(x_isolated, i).get(dim) -= gsl::at(coords_isolated, i)[dim];
484 : }
485 : gsl::at(euclidean_distance, i) = get(magnitude(gsl::at(x_isolated, i)));
486 : if (falloff_widths_.has_value()) {
487 : gsl::at(windows, i) = exp(-square(gsl::at(euclidean_distance, i)) /
488 : square(gsl::at(*falloff_widths_, i)));
489 : } else {
490 : gsl::at(windows, i) = make_with_value<DataVector>(x, 1.);
491 : }
492 : gsl::at(isolated_vars, i) = get_isolated_vars<requested_superposed_tags>(
493 : *gsl::at(superposed_objects_, i), gsl::at(x_isolated, i));
494 : }
495 : auto flat_vars = flatness_.variables(x, requested_superposed_tags{});
496 : typename VarsComputer::Cache cache{get_size(*x.begin())};
497 : const VarsComputer computer{std::move(mesh),
498 : std::move(inv_jacobian),
499 : x,
500 : angular_velocity_,
501 : expansion_,
502 : linear_velocity_,
503 : falloff_widths_,
504 : std::move(x_isolated),
505 : std::move(windows),
506 : std::move(flat_vars),
507 : std::move(isolated_vars)};
508 : return {cache.get_var(computer, RequestedTags{})...};
509 : }
510 :
511 : template <typename TagsList, typename... Args>
512 0 : tuples::tagged_tuple_from_typelist<TagsList> get_isolated_vars(
513 : const IsolatedObjectBase& isolated_object, const Args&... args) const {
514 : return call_with_dynamic_type<tuples::tagged_tuple_from_typelist<TagsList>,
515 : IsolatedObjectClasses>(
516 : &isolated_object, [&args...](const auto* const derived) {
517 : return derived->variables(args..., TagsList{});
518 : });
519 : }
520 : };
521 :
522 : /// \cond
523 : template <typename IsolatedObjectBase, typename IsolatedObjectClasses>
524 : PUP::able::PUP_ID Binary<IsolatedObjectBase, IsolatedObjectClasses>::my_PUP_ID =
525 : 0; // NOLINT
526 : /// \endcond
527 :
528 : } // namespace Xcts::AnalyticData
|
