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             : #include "NumericalAlgorithms/SphericalHarmonics/Strahlkorper.hpp"
      13             : 
      14           1 : namespace evolution::Ringdown {
      15             : 
      16             : /*!
      17             :  * \brief This function finds a safe ringdown excision radius for starting
      18             :  * the ringdown of a common horizon from a binary inspiral. It does this by
      19             :  * finding a radius that will enclose every point on strahlkorpers that
      20             :  * represent excisions A/B from the inspiral at the match time. It does this by
      21             :  * taking excisions A/B from the inertial frame and transforming them to the
      22             :  * ringdown grid frame. It iterates over all the points from excisions A/B in
      23             :  * the ringdown-grid-frame until it finds the point that's the farthest from the
      24             :  * geometric center of the ringdown excision. This is now the minimum radius we
      25             :  * can choose to enclose the excisions A/B from the inspiral. This is not the
      26             :  * ideal radius however, we then choose a radius that's 3/4 of the way between
      27             :  * the radius of the common horizon at the match time and the minimum radius
      28             :  * that encloses excisions A/B from the inspiral and check for multiple l_max
      29             :  * values of excisions A/B to ensure they fit and the ringdown can start.
      30             : 
      31             :  * \details It does this by constructing inspiral-grid-frame AhA/AhB excision
      32             :  * strahlkorpers from the radii and centers passed to this function. It then
      33             :  * maps those excisions to the inertial-frame using the inspiral domain
      34             :  * and functions of time from the inspiral volume data and subfile supplied. We
      35             :  * then construct a test ringdown domain that has all the corrected functions of
      36             :  * time from ComputeRingdownShapeAndTranslationFoT.py and an initial guess for
      37             :  * the inner radius. More details are outlined below.
      38             :  * \param path_to_volume_data The full path to the volume data containing the
      39             :  * domain and functions of time
      40             :  * \param volume_subfile_name Subfile containing volume data output from the
      41             :  * inspiral
      42             :  * \param path_to_horizons_h5 Path to h5 file containing horizon data for AhA/B
      43             :  * \param surface_subfile_name Subfile containing horizon data for AhA/B
      44             :  * \param path_to_AhC_distorted_h5 Path to h5 file containing ringdown shape
      45             :  * coefficients computed using ComputeRingdownShapeAndTranslationFoT.py
      46             :  * \param AhC_distorted_subfile_names Subfiles in the h5 file containing
      47             :  * shape coefficients
      48             :  * \param match_time The time to match the functions of time
      49             :  * \param settling_timescale Timescale at which the functions of time settle to
      50             :  * constant values
      51             :  * \param excision_A_radius The radius of excision A from the inspiral grid
      52             :  * frame
      53             :  * \param excision_B_radius The radius of excision B from the inspiral grid
      54             :  * frame
      55             :  * \param excision_A_center The center of excision A from the inspiral grid
      56             :  * frame
      57             :  * \param excision_B_center The center of excision B from the inspiral grid
      58             :  * frame
      59             :  * \param exp_func_and_2_derivs Expansion function of time from the inspiral
      60             :  * \param exp_outer_bdry_func_and_2_derivs Outer boundary expansion function of
      61             :  * time from the inspiral
      62             :  * \param rot_func_and_2_derivs Rotation function of time from the inspiral
      63             :  * \param trans_func_and_2_derivs The corrected translation function of time
      64             :  * computed by ComputeRingdownShapeAndTranslationFoT.py
      65             :  * \param match_time_tol The difference allowed between the match time requested
      66             :  * and the time found in h5 files
      67             :  *
      68             :  * Using the corrected functions of time and our initial guess for the ringdown
      69             :  * inner radius (the excision radius), excisions A/B are transformed from the
      70             :  * inertial frame to the ringdown grid frame. The inner radius is iterated upon
      71             :  * using 2 main loops, the outer loop which changes the L_max of the excisions
      72             :  * A/B being transformed from the inspiral and the inner loop that changes the
      73             :  * excision radius used in the ringdown. The general flow for these loops is as
      74             :  * follows:
      75             :  *
      76             :  * Let $ahc\_average\_radius$ be the average radius of the common horizon in the
      77             :  * inertial frame and since have strahlkorper data from the inertial frame AhC
      78             :  * we set $ahc\_average\_radius$ so that lambda00=0 at the match time.
      79             :  *
      80             :  * Let $ringdown\_excision\_radius$ be the radius of (spherical) excision
      81             :  * boundary in ringdown grid frame.
      82             :  *
      83             :  * Let $ringdown\_excision\_factor$ be the factor we multiply the
      84             :  * $ahc\_average\_radius$ by to get the $ringdown\_excision\_radius$.
      85             :  *
      86             :  * The goal is to find a good $ringdown\_excision\_factor$
      87             :  *
      88             :  * 1) Choose an initial guess for the $ringdown\_excision\_factor$. We choose
      89             :  * this to be 0.94 by default.
      90             :  *
      91             :  * 2) Choose an initial guess for current_l_max, the angular resolution of
      92             :  * spherical Strahlkorpers that match excision A and B (not horizons) in the
      93             :  * inspiral grid frame. We choose current_l_max=20 intially.
      94             :  *
      95             :  * 3) Construct Strahlkorpers in the Inspiral grid frame that correspond to
      96             :  * excisions A/B, the excision regions (not the horizons) of the two individual
      97             :  * holes. These are spherical in the inspiral grid frame, and have
      98             :  * L=current_l_max.
      99             :  *
     100             :  * 4) Map the Strahlkorper excisions A/B to the inertial frame. Now they are not
     101             :  * spherical.
     102             :  *
     103             :  * 4a) Using the current $ringdown\_excision\_factor$, construct a test domain
     104             :  * to map excisions A/B from the inertial frame to the ringdown grid frame. This
     105             :  * domain's functions of time should contain the scaling and rotation functions
     106             :  * of time from the inspiral that settle to const and the shape and corrected
     107             :  * translation function of time output by
     108             :  * ComputeRingdownShapeAndTranslationFoT.py
     109             :  *
     110             :  * 4b) Map the Strahlkorper excisions A/B to the Ringdown grid frame. They are
     111             :  * still not spherical.
     112             :  *
     113             :  * 4c) Loop through all the points on Ringdown grid excisions A/B, and find the
     114             :  * point that is the maximum distance from the origin. Call that maximum
     115             :  * distance $min\_excision\_radius$.
     116             :  *
     117             :  * 4d) A Ringdown grid sphere of radius $min\_excision\_radius$ centered at the
     118             :  * origin will (barely) enclose both excisions A/B, so choose
     119             :  * $minimum\_ringdown\_excision\_factor = \frac{min\_excision\_radius}
     120             :  * {ahc\_average\_radius}$ as the minimum possible $ringdown\_excision\_factor$.
     121             :  * And then choose the new value of $ringdown\_excision\_factor$ to be
     122             :  * \f{equation}{
     123             :  * ringdown\_excision\_factor=1-0.25*(1-minimum\_ringdown\_excision\_factor)
     124             :  * \f}
     125             :  * Which puts the excision radius 3/4 of the way between the
     126             :  * $ahc\_average\_radius$ and the $min\_excision\_radius$
     127             :  *
     128             :  * 4e) If this is the first inner iteration, or if
     129             :  * $|ringdown\_excision\_factor-previous\_ringdown\_excision\_factor|>eps$, then
     130             :  * set $previous\_ringdown\_excision\_factor=ringdown\_excision\_factor$ and
     131             :  * goto 3a, where 3a is repeated using the new $ringdown\_excision\_factor$. We
     132             :  * choose $eps = \frac{10^{-3}}{q}$ where q is an approximation of the mass
     133             :  * ratio. The reason this is iterative is because the new
     134             :  * $ringdown\_excision\_factor$ in 3d depends on the ringdown_excision_factor
     135             :  * used when constructing the test domain.
     136             :  *
     137             :  * 5) If we get here, then inner iteration has converged.  But that iteration
     138             :  * depends on the $current\_l\_max$ used to construct excisions A/B in step 2.
     139             :  * If this is the first outer iteration, or if
     140             :  * $|ringdown\_excision\_factor-previous\_ringdown\_excision\_factor|>eps$, then
     141             :  * set $previous\_ringdown\_excision\_factor=ringdown\_excision\_factor$, and
     142             :  * increment $current\_l\_max$ by some increment (which we choose to be 6), and
     143             :  * go back to 2.
     144             :  *
     145             :  * 6) If we get here, the outer iteration is converged and now we have a final
     146             :  * $ringdown\_excision\_factor$.
     147             :  *
     148             :  * \note This implementation does not do everything SpEC does yet, it is
     149             :  * currently missing rescaling the shape coefficients held in
     150             :  * 'path_to_AhC_distorted_h5' by $excision\_radius / average\_ahc\_radius$ every
     151             :  * time it constructs a test ringdown domain. This step helps the shape of the
     152             :  * excision match the shape of the apparent horizon.
     153             :  */
     154           1 : double minimum_ahc_excision_radius(
     155             :     const std::string& path_to_volume_data,
     156             :     const std::string& volume_subfile_name,
     157             :     const std::string& path_to_horizons_h5,
     158             :     const std::string& surface_subfile_name,
     159             :     const std::string& path_to_AhC_distorted_h5,
     160             :     const std::vector<std::string>& AhC_distorted_subfile_names,
     161             :     double match_time, double settling_timescale, double excision_A_radius,
     162             :     double excision_B_radius, std::array<double, 3> excision_A_center,
     163             :     std::array<double, 3> excision_B_center,
     164             :     const std::optional<std::array<double, 3>>& exp_func_and_2_derivs =
     165             :         std::nullopt,
     166             :     const std::optional<std::array<double, 3>>&
     167             :         exp_outer_bdry_func_and_2_derivs = std::nullopt,
     168             :     const std::optional<std::vector<std::array<double, 4>>>&
     169             :         rot_func_and_2_derivs = std::nullopt,
     170             :     const std::optional<std::array<std::array<double, 3>, 3>>&
     171             :         trans_func_and_2_derivs = std::nullopt,
     172             :     double match_time_tol = 1e-12);
     173             : }  // namespace evolution::Ringdown