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/Strahlkorper/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
|
