SpECTRE Documentation Coverage Report
Current view: top level - PointwiseFunctions/AnalyticData/Xcts - Binary.hpp Hit Total Coverage
Commit: 047f5140268fc2680430d095fc7d96883747c68e Lines: 4 59 6.8 %
Date: 2024-12-09 22:54:54
Legend: Lines: hit not hit

          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

Generated by: LCOV version 1.14