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 <cstddef>
       8             : #include <string>
       9             : #include <vector>
      10             : 
      11             : #include "DataStructures/DataVector.hpp"
      12             : 
      13             : /*!
      14             :  * \brief Functionality for evolving a ringdown following a compact-binary
      15             :  * merger.
      16             :  */
      17           1 : namespace evolution::Ringdown {
      18             : /*!
      19             :  * \brief This function is used to transition from inspiral to ringdown. It
      20             :  * reads inertial frame common horizon Strahlkorper coefs from a file and
      21             :  * returns the Strahlkorper's temporary frame coefs and the inertial frame
      22             :  * geometric centers at multiple times which will be used to initialize the
      23             :  * shape and translation function of time for the ringdown.
      24             :  *
      25             :  * \details Reads common horizon Strahlkorpers (assumed to be in the inertial
      26             :  * frame) from a file, then transforms them to a temporary ringdown domain
      27             :  * defined by the expansion, rotation, and translation maps from the inspiral
      28             :  * specified by `exp_func_and_2_derivs`, `exp_outer_bdry_func_and_2_derivs`,
      29             :  * `rot_func_and_2_derivs`, and `trans_func_and_2_derivs`. The expansion and
      30             :  * rotation functions of time from the inspiral are the same as the ringdown
      31             :  * frame's expansion and rotation maps at the given `match_time`, but later
      32             :  * settle to constant values by the given `settling_timescale`. The translation
      33             :  * function of time supplied from the inspiral is not the translation map we'll
      34             :  * use in the final ringdown, it is only used to correctly map the common
      35             :  * horizon's geometric center so that we can make the corrected translation map
      36             :  * for the final ringdown. How to correct the translation map is explained in
      37             :  * detail below. A shape map is not specified, because we do not yet know the
      38             :  * shape coefficients of the common horizon for the ringdown. There are 4 frames
      39             :  * we transition between to find the shape coefficients and correct the
      40             :  * translation map. The frames are as follows:
      41             :  *
      42             :  * $\textbf{Inertial frame}$: This is the same for ringdown and inspiral, by
      43             :  * definition. It's the frame where we measure the waveforms, and it's the frame
      44             :  * in which we define the tensor components (e.g. the 'x' and 'y' in 'g_xy').
      45             :  *
      46             :  * $\textbf{Ringdown Grid Frame}$: This is the frame where the ringdown excision
      47             :  * boundary is spherical and where the ringdown grid lives.
      48             :  *
      49             :  * $\textbf{Temporary frame}$: The inertial frame pushed through the inverse
      50             :  * translation-rotation-expansion map, where "translation" is the inspiral's
      51             :  * translation map, and "rotation" and "expansion" are new rotation and
      52             :  * expansion maps that match the inspiral's rotation and expansion at the outer
      53             :  * boundary at the match time, have some smooth falloff in the interior, and
      54             :  * have a "settle to constant" time dependence. This frame is almost the
      55             :  * ringdown distorted frame, except for an uncorrected translation map which is
      56             :  * explained in detail below.
      57             :  *
      58             :  * $\textbf{Ringdown distorted frame}$: This is the ringdown grid frame pushed
      59             :  * through the (final) forward distortion map. Or equivalently, this is the
      60             :  * inertial frame pushed through the inverse translation-rotation-expansion map,
      61             :  * where here "translation" means the final corrected translation map, not the
      62             :  * inspiral's translation map, and "rotation" and "expansion" are the same as in
      63             :  * the temporary frame.
      64             :  *
      65             :  * We get the shape coefficients by transforming the common horizon to the
      66             :  * temporary frame which is almost the ringdown distorted frame except for the
      67             :  * uncorrected translation map. The translation map is corrected by transforming
      68             :  * the common horizon from the temporary frame to the inertial frame and saving
      69             :  * the geometric center points at multiple times. We then take those geometric
      70             :  * center points and build the corrected translation function of time. This is
      71             :  * done because the center of the Strahlkorper in this temporary frame is NOT
      72             :  * the origin (this is the case whether or not you say Recenter=true when
      73             :  * transforming the Strahlkorper), but the distortion map does its distortion
      74             :  * about the origin.
      75             :  * Possible ways to account for this:
      76             :  *
      77             :  * 1) Put a translation map before the distortion map.
      78             :  *
      79             :  * 2) Change the center of the distortion map.  (then we need to live with
      80             :  * this during the ringdown).
      81             :  *
      82             :  * 3) Correct the current translation map so that the excision boundary maps to
      83             :  * the correct place.
      84             :  *
      85             :  * We choose 3). Section 6 of \cite Hemberger2012jz explains in more
      86             :  * detail of how we initialize the shape/translation map, but idea is that in
      87             :  * Eq. (104), the horizon can be written as
      88             :  * \f{equation}{
      89             :  * x^{\bar{i}}_{AH} = x^{\bar{i}}_{AHc} + n^{\bar{i}} \sum_{lm}S_{lm}
      90             :  * Y_{lm}(\theta,\phi)
      91             :  * \f}
      92             :  * where $x^{\bar{i}}_{AH}$ is the common horizon Strahklorper in the temporary
      93             :  * frame, $x^{\bar{i}}_{AHc}$ is the expansion center of the Strahlkorper in the
      94             :  * temporary frame, which is time-dependent and not zero, $n^{\bar{i}}$ is the
      95             :  * direction unit vector in the ($\theta,\phi$) direction relative to
      96             :  * $x^{\bar{i}}_{AHc}$, $\bar{i}$ is the index corresponding to the
      97             :  * temporary frame, and $S_{lm}$ are the Strahlkorper coefficients.
      98             :  * Now,
      99             :  * \f{equation}{
     100             :  * x^i = T0^i + M^i_{\bar{i}} x^{\bar{i}}
     101             :  * \f}
     102             :  * where $T0^i$ is the inspiral translation map, and $M^i_{\bar{i}}$ is
     103             :  * scaling+rotation. Thus,
     104             :  * \f{equation}{
     105             :  * x^i_AH = T0^i + M^i_{\bar{i}} x^{\bar{i}}_{AH}
     106             :  *        = T0^i + M^i_{\bar{i}} x^{\bar{i}}_{AHc} + M^i_{\bar{i}} n^{\bar{i}}
     107             :  *          \sum_{lm} S_{lm} Y_{lm}
     108             :  * \f}
     109             :  * Therefore if you define a new translation map $T^i$ as in Eq. (107)
     110             :  * (corrected translation map that will be used in the final ringdown)
     111             :  * $T^i = T0^i + M^i_{\bar{i}} x^{\bar{i}}_{AHc}$ (that is, you define $T^i$ to
     112             :  * be the same as $x^i_{AHc}$), then you can rewrite the relationship as
     113             :  * \f{equation}{
     114             :  * x^{i}_{AH} = T^i + M^i_{\bar{i}} \sum_{lm} S_{lm} Y_{lm} n^{\bar{i}}
     115             :  * \f}
     116             :  * Therefore we use a new map $x^i = T^i + M^i_{\bar{i}} x^{\tilde{i}}$ where
     117             :  * $\tilde{i}$ refers to the final ringdown distorted frame, $x^{\tilde{i}}$ is
     118             :  * a new coordinate where $x^{\tilde{i}}_{AHc} = 0$ and the coefficients of the
     119             :  * AH in the $x^{\tilde{i}}$ frame can be used unchanged (except for a minus
     120             :  * sign) in the distortion map that connects $x^{i_{grid}}$ and $x^{\tilde{i}}$.
     121             :  * This new map for $x^i$ ensures that the excision boundary is mapped to the
     122             :  * correct position.
     123             :  * \note Only temporary frame common horizon coefs and inertial frame geometric
     124             :  * center points within `requested_number_of_times_from_end` times from the
     125             :  * final time are returned. The functions of time are computed by
     126             :  * ComputeRingdownShapeAndTranslationFoT.py using the output from this function
     127             :  */
     128             : std::pair<std::vector<DataVector>, std::vector<std::array<double, 3>>>
     129           1 : strahlkorper_coefs_and_centers(
     130             :     const std::string& path_to_volume_data,
     131             :     const std::string& volume_subfile_name,
     132             :     const std::string& path_to_horizons_h5,
     133             :     const std::string& surface_subfile_name,
     134             :     size_t requested_number_of_times_from_end, double match_time,
     135             :     double settling_timescale,
     136             :     const std::optional<std::array<double, 3>>& exp_func_and_2_derivs =
     137             :         std::nullopt,
     138             :     const std::optional<std::array<double, 3>>&
     139             :         exp_outer_bdry_func_and_2_derivs = std::nullopt,
     140             :     const std::optional<std::vector<std::array<double, 4>>>&
     141             :         rot_func_and_2_derivs = std::nullopt,
     142             :     const std::optional<std::array<std::array<double, 3>, 3>>&
     143             :         trans_func_and_2_derivs = std::nullopt);
     144             : 
     145             : // Hardcoded value used in SpEC
     146           0 : constexpr double expansion_center_tolerance = 1e-8;
     147             : }  // namespace evolution::Ringdown