diff --git a/docs/tutorial_abd.rst b/docs/tutorial_abd.rst index 8136a26a..2fa76ebe 100644 --- a/docs/tutorial_abd.rst +++ b/docs/tutorial_abd.rst @@ -74,7 +74,7 @@ relativists. First, the ABD class uses the Moreschi-Boyle convention, which is * the standard convention used in NR. See Appendices B and C of `arXiv:2010.15200 `_ for more information. Second, the Newman-Penrose shear :math:`\sigma` is *not* the gravitational-wave strain :math:`h`. -In the Moreschi-Boyle convention, we have the relation :math:`h = \bar{\sigma}`, but +In the Moreschi-Boyle convention, we have the relation :math:`h = 2\bar{\sigma}`, but this relation is not gauranteed in every convention and only holds for asymptotic quantities. @@ -277,3 +277,81 @@ These components of a BMS transformation can also all be stored in the :meth:`scri.bms_transformations.BMSTransformation` class, which allows for things like reording the components of the BMS transformation, inverting BMS transformations, and composing BMS transformations. For more, see :ref:`bms_transformations`. + +========== +BMS Frames +========== + +All waveforms at future null infinity (and all waveforms more generally) are functions of coordinates. +Therefore, there are certain "frames" which may be more useful than others, like that of a rest frame. +For waveforms at future null infinity, the number of coordinate freedoms, i.e., the symmetries, that they +exhibit is infinite and is summarized by a group known as the BMS group. This controls the types of frames +that one may map waveforms to. Because GR is covariant, there is no preferred frame. However, for performing +analysis on waveforms or building waveform models, it turns out that there are certain frames that are +more useful than others. In particular, within GR one can extend the notion of a rest frame to something called +a "superrest frame" (see arXiv:2405.08868 or arXiv:2208.04356 for more details), which typically yields waveforms +that are easier to understand/analyze. Effectively, mapping to this frame amounts to mapping the system to be +in the center-of-mass frame, with no instananeous memory, and its angular velocity in the z-direction. For example, +for a remnant black hole, this corresponds to making the coordinates match those of the usual Kerr metric and is +therefore incredibly useful (and necessary) for fitting QNMs to NR waveforms. + +The function ``scri.asymptotic_bondi_data.map_to_superrest_frame`` maps to this exact frame. +In particular, it takes as input: + +* ``t_0``, the time at which to map to the superrest frame; + +* ``target_PsiM_input``, the target Moreschi supermomentum; this should be ``None`` to map to the superrest frame, + but to map to the PN BMS frame one should input the PN Moreschi supermomentum (see arXiv:2208.04356). + +* ``target_strain_input``, the target strain; this should be ``None`` to map to the superrest frame, + but to map to the PN BMS frame one should input the PN strain (see arXiv:2208.04356). + +* ``padding_time``, the time window about ``t_0`` to be used when finding the BMS transformation to the superrest frame. + +======================================================== +Creating an AsymptoticBondiData Object from a CCE Output +======================================================== + +For processing the output of SpECTRE CCE, one may use the function ``scri.SpEC.file_io.create_abd_from_h5``. +This function takes as input the path to SpECTRE CCE's output file (via the option ``file_name``) and +creates an ``abd`` object from said file. +It also can perform a number of other important post-processing steps, such as: + +* time translate the time array of the waveforms by the radius of the worldtube; this ensures that the CCE waveforms + are more closely aligned (in time) with extrapolated waveform. This is performed via the ``radius`` option. + +* scale out the total Christoudoulou mass of the system from each waveform. This is performed via the ``ch_mass`` option. + +* interpolate to a coarser time array, such as the time array of the worldtube. This is performed via the ``t_interpolate`` option. + +* map to the superrest BMS frame at some time. This is performed via the ``t_0_superrest`` and ``padding_time`` options. + E.g., to make reasonable-looking waveforms, one should map to the superrest frame at some time after junk radiation; + ``t_0_superrest`` is the time at which to map to this frame, and ``padding_time`` is the window around ``t_0_superrest`` + that is used when computing this BMS transformation. ``t_0_superrest - padding_time`` should be after junk radiation. + ``padding_time`` should be a few hundred :math:`h=2\overline{\sigma}` (where :math:`\sigma` is the shear), e.g., two orbital periods. + The function used to do this is ``abd.map_to_superrest_frame`` (see the "BMS Frames" section). + +We recommend including all of these post-processing steps when processing SpECTRE CCE output. + +To obtain the strain :math:`h` from the ``abd`` object, one can use the function ``scri.asymptotic_bondi_data.map_to_superrest_frame.MT_to_WM`` via ``h = MT_to_WM(2.0*abd.sigma.bar)``. +This is because the strain :math:`h` is related to the shear :math:`\sigma` via :math:`h=2\overline{\sigma}`. + +Example usage of this function could be: + +.. code-block:: python + + >>> import scri + >>> import matplotlib.pyplot as plt + >>> abd = scri.SpEC.file_io.create_abd_from_h5( + file_name='CharacteristicExtractVolume_R0292.h5', + file_format="spectrecce", + ch_mass=1.0, + t_interpolate=t_worldtube, + t_0_superrest=1600, + padding_time=200 + ) + >>> h = abd.h + >>> plt.plot(h.t, h.data[:, h.index(2,2)]) + >>> plt.show() + +For more on the :meth:`scri.WaveformModes` class, i.e., what :math:`h` is in the above code, see https://github.com/moble/scri/blob/main/docs/tutorial_waveformmodes.rst. diff --git a/scri/SpEC/file_io/__init__.py b/scri/SpEC/file_io/__init__.py index 4eb185cb..a1cff5aa 100644 --- a/scri/SpEC/file_io/__init__.py +++ b/scri/SpEC/file_io/__init__.py @@ -33,7 +33,6 @@ rpxmb = rotating_paired_xor_multishuffle_bzip2 - def translate_data_types_GWFrames_to_waveforms(d): if d < 8: return {0: UnknownDataType, 1: h, 2: hdot, 3: psi4, 4: psi3, 5: psi2, 6: psi1, 7: psi0}[d] @@ -149,7 +148,6 @@ def read_from_h5(file_name, **kwargs): pass # Probably couldn't find the metadata.json/metadata.txt file try: # Make sure the h5py.File gets closed, even in the event of an exception - # Add the old history to the new history, if found try: try: @@ -420,7 +418,48 @@ def write_to_h5(w, file_name, file_write_mode="w", attributes={}, use_NRAR_forma f.close() -def create_abd_from_h5(file_format, convention="SpEC", **kwargs): +def make_variable_dimensionless(WM, ch_mass=1.0): + """ + Make variable dimensionless by scaling out Christodoulou mass. + + Parameters: + ----------- + WM : scri.WaveformModes + Waveform data. + ch_mass : float, optional + Total Christodoulou mass of the system. + [Default: 1.0]. + + """ + + if WM.m_is_scaled_out: + print("Data is already dimensionless!") + return + + if WM.dataType in [psi4, psi3, psi2, psi1, psi0]: + unit_scale_factor = (ch_mass) ** (WM.dataType - 4) + elif WM.dataType == h: + unit_scale_factor = 1 / ch_mass + elif WM.dataType == hdot: + unit_scale_factor = 1.0 + else: + raise ValueError("DataType not determined.") + + WM.t = WM.t / ch_mass + WM.data = WM.data * unit_scale_factor + WM.m_is_scaled_out = True + + +def create_abd_from_h5( + file_format, + convention="SpEC", + radius=None, + ch_mass=None, + t_interpolate=None, + t_0_superrest=None, + padding_time=None, + **kwargs, +): """Returns an AsymptoticBondiData object with waveform data loaded from specified H5 files. The AsymptoticBondiData class internally uses the Moreschi-Boyle conventions, see the following reference: @@ -429,18 +468,28 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): Parameters ---------- - file_format : 'SXS', 'CCE', 'SpECTRECCE1', 'RPDMB', or 'RPXMB' + file_format : 'SXS', 'SpECTRECCE_v1', 'RPDMB', or 'RPXMB' The H5 files may be in the one of the following file formats: * 'SXS' - Dimensionless extrapolated waveform files found in the SXS Catalog, also known as NRAR format. - * 'CCE' - Asymptotic waveforms output by SpECTRE CCE. These are not dimensionless. - * 'SpECTRECCE1' - Asymptotic waveforms output by SpECTRE CCE's original format. - * SpECTRECCE' - Asymptotic waveforms output by SpECTRE CCE's newest format. (NOTE: this may break compatibility) + * 'SpECTRECCE_v1' - Asymptotic waveforms output by SpECTRE CCE v1. * 'RPDMB' - Dimensionless waveforms compressed using the rotating_paired_diff_multishuffle_bzip2 format. * 'RPXMB' - Dimensionless waveforms compressed using the rotating_paired_xor_multishuffle_bzip2 format. convention : 'SpEC' or 'Moreschi-Boyle' The data conventions of the waveform data that will be loaded. This defaults to 'SpEC' since this will be - most often used with 'SpEC' convention waveforms. - + most often used with 'SpEC' convention waveforms. The output convention is Moreschi-Boyle. + radius : str, optional + Worldtube radius used when running CCE; only needed for versions of SpECTRE before PR #5985. + The time array of the worldtube is translated by this radius. + ch_mass : float, optional + Total Christodoulou mass of the system. + t_interpolate : float array, optional + Time array to interpolate to, e.g., the time array of the worldtube. + t_0_superrest : float, optional + When to map to the BMS superrest frame. + Typically a few hundred M after the junk radiation is sufficient. + padding_time : float, optional + Time window length around t_0_superrest to use when mapping to the superrest frame. + Typically a few hundred M or a few orbits is sufficient. Keyword Parameters ------------------ @@ -466,14 +515,28 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): # Load waveform data from H5 files into WaveformModes objects WMs = {} filenames = {} - if file_format == "spectrecce" or file_format == "spectrecce1": - file_name = kwargs.pop("file_name") + if file_format == "spectrecce_v1": + try: + file_name = kwargs.pop("file_name") + except: + raise ValueError('Need to specify "file_name" option!') + with h5py.File(file_name, "r") as f: - cce = f["Cce"] - time = cce["Strain.dat"][:, 0] + for x in f.keys(): + if "Spectre" in x: + cce_key = x + radius = cce_key.split("R")[1][:4] + data_label_suffix = "" + break + else: + cce_key = "Cce" + data_label_suffix = ".dat" + + cce = f[cce_key] + time = cce[f"Strain{data_label_suffix}"][:, 0] indices = monotonic_indices(time) time = time[indices] - ell_max = int(np.sqrt((cce["Strain.dat"].shape[1]-1) / 2) - 1) + ell_max = int(np.sqrt((cce[f"Strain{data_label_suffix}"].shape[1] - 1) / 2) - 1) for data_label in ["Psi4", "Psi3", "Psi2", "Psi1", "Psi0", "Strain"]: if data_label != "Strain": dataType = DataNames.index(data_label) @@ -481,7 +544,7 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): dataType = DataNames.index("h") WMs[data_label] = WaveformModes( t=time, - data=cce[f"{data_label}.dat"][indices, 1:].view(np.complex128), + data=cce[f"{data_label}{data_label_suffix}"][indices, 1:].view(np.complex128), ell_min=0, ell_max=ell_max, frameType=Inertial, @@ -508,13 +571,11 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): WMs[data_label] = rotating_paired_xor_multishuffle_bzip2.load(filenames[data_label])[0] WMs[data_label].to_inertial_frame() elif file_format == "rpdmb" or file_format == "rpdm": - WMs[data_label] = WaveformModes.from_sxs( - sxs.rpdmb.load(filenames[data_label]) - ) + WMs[data_label] = WaveformModes.from_sxs(sxs.rpdmb.load(filenames[data_label])) else: raise ValueError( f"File format '{file_format}' not recognized. " - "Must be 'SXS', 'CCE', 'SpECTRECCE', 'RPDMB', or 'RPXMB'." + "Must be 'SXS', 'CCE', 'SpECTRECCE_v1', 'RPDMB', or 'RPXMB'." ) if kwargs: @@ -525,6 +586,7 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): # Sanity check if not WMs: raise ValueError("No filenames have been provided. The data of at least one waveform quantity is required.") + WM_ref = WMs[list(WMs.keys())[0]] for i in WMs: if not (WM_ref.t == WMs[i].t).all(): @@ -533,6 +595,15 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): f"for {list(WMs.keys())[i].data_type_string} has a different set of times." ) + for i in WMs: + # Make waveforms dimensionless (if they already are, does nothing) + if ch_mass is not None: + make_variable_dimensionless(WMs[i], ch_mass) + + # Adjust time by worldtube radius + if file_format == "spectrecce_v1": + WMs[i].t -= float(radius) + # Create an instance of AsymptoticBondiData abd = AsymptoticBondiData( time=WM_ref.t, @@ -581,4 +652,14 @@ def create_abd_from_h5(file_format, convention="SpEC", **kwargs): ) abd.sigma = abd.sigma.bar + # Interpolate to finer time array, if specified + if t_interpolate is not None: + idx1 = np.argmin(abs(t_interpolate - abd.t[0])) + 1 + idx2 = np.argmin(abs(t_interpolate - abd.t[-1])) + 1 - 1 + abd = abd.interpolate(t_interpolate[idx1:idx2]) + + # Map to superrest frame at some time over some window, if specified + if t_0_superrest is not None and padding_time is not None: + abd, BMS, _ = abd.map_to_superrest_frame(t_0=t_0_superrest, padding_time=padding_time) + return abd diff --git a/scri/__init__.py b/scri/__init__.py index 2551d236..9cb992d8 100644 --- a/scri/__init__.py +++ b/scri/__init__.py @@ -159,6 +159,7 @@ def version_info(): from . import sample_waveforms, SpEC, LVC, utilities from .SpEC import rpxmb +from .SpEC.file_io import create_abd_from_h5 from .bms_transformations import LorentzTransformation, BMSTransformation diff --git a/scri/asymptotic_bondi_data/__init__.py b/scri/asymptotic_bondi_data/__init__.py index ac3aaa39..da617755 100644 --- a/scri/asymptotic_bondi_data/__init__.py +++ b/scri/asymptotic_bondi_data/__init__.py @@ -1,7 +1,9 @@ import numpy as np from spherical_functions import LM_total_size +from .. import WaveformModes from .. import ModesTimeSeries from .. import Inertial +from .. import h as h_DataType class AsymptoticBondiData: @@ -114,6 +116,20 @@ def sigma(self, sigmaprm): self._sigma[:] = sigmaprm return self.sigma + @property + def h(self): + h_mts = 2.0 * self._sigma.bar + return WaveformModes( + t=h_mts.t, + data=np.array(h_mts)[:, h_mts.index(abs(h_mts.s), -abs(h_mts.s)) :], + ell_min=abs(h_mts.s), + ell_max=h_mts.ell_max, + frameType=Inertial, + dataType=h_DataType, + r_is_scaled_out=True, + m_is_scaled_out=True, + ) + @property def psi4(self): return self._psi4