Working with the Public ICCUB Catalog and Waveform Integration¶

This tutorial demonstrates how to work with waveforms from the public ICCUB (Institute of Cosmos Sciences) numerical relativity catalog, hosted at egrav.icc.ub.edu, and how to perform waveform integration from Ψ₄ (Weyl scalar) to strain h.

ICCUB simulations provide Ψ₄ data, which needs to be integrated twice to obtain the gravitational wave strain. PyART provides tools for this integration, and also stores a pre-computed strain in the downloaded file that we can use as a reference.

Setup¶

%matplotlib inline
%config InlineBackend.figure_format = 'retina'

import matplotlib.pyplot as plt
import numpy as np
from PyART.catalogs.icc_public import Waveform_ICC
from PyART.analysis.integrate_wave import IntegrateMultipole

Load ICCUB Waveform¶

Waveform_ICC downloads the simulation from the public catalog the first time it’s requested (download=True), and reuses the cached files under path afterwards. ellmax is capped at the highest multipole the public release provides (l<=4).

# Public ICCUB catalog entry to use
path = './local_data/icc/'
ID = 1

icc = Waveform_ICC(path=path, ID=ID, download=True, ellmax=4)
icc.load_psi4lm()  # psi4 modes are not loaded automatically

print("Waveform loaded successfully!")

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[2], line 5
# Public ICCUB catalog entry to use
path = './local_data/icc/'
ID = 1

----> 5 icc = Waveform_ICC(path=path, ID=ID, download=True, ellmax=4)
icc.load_psi4lm()  # psi4 modes are not loaded automatically

print("Waveform loaded successfully!")

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/PyART/catalogs/icc_public.py:78, in Waveform_ICC.__init__(self, path, ID, download, load, ellmax, nu_rescale, extraction)
json_catalog = "egrav_data.json"
if not os.path.exists(json_catalog):
---> 78     fetch_and_save_egrav_json()
self.ID = ID
self.path = path

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/PyART/catalogs/icc_public.py:26, in fetch_and_save_egrav_json(output_path)
"""
Download the JSON catalog of public ICC waveforms and save it locally.
Parameters
   (...)     23     List of entries in the catalog.
"""
url = "https://egrav.icc.ub.edu/site/list-grid"
---> 26 response = requests.get(url, verify=False)
response.raise_for_status()
data = response.json()

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/requests/api.py:87, in get(url, params, **kwargs)
def get(
   url: _t.UriType, params: _t.ParamsType = None, **kwargs: Unpack[_t.GetKwargs]
) -> Response:
   r"""Sends a GET request.

   :param url: URL for the new :class:`Request` object.
   (...)     84     :rtype: requests.Response
   """
---> 87     return request("get", url, params=params, **kwargs)

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/requests/api.py:71, in request(method, url, **kwargs)
# By using the 'with' statement we are sure the session is closed, thus we
# avoid leaving sockets open which can trigger a ResourceWarning in some
# cases, and look like a memory leak in others.
with sessions.Session() as session:
---> 71     return session.request(method=method, url=url, **kwargs)

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/requests/sessions.py:651, in Session.request(self, method, url, params, data, headers, cookies, files, auth, timeout, allow_redirects, proxies, hooks, stream, verify, cert, json)
send_kwargs = {
   "timeout": timeout,
   "allow_redirects": allow_redirects,
}
send_kwargs.update(settings)
--> 651 resp = self.send(prep, **send_kwargs)
return resp

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/requests/sessions.py:784, in Session.send(self, request, **kwargs)
start = preferred_clock()
# Send the request
--> 784 r = adapter.send(request, **kwargs)
# Total elapsed time of the request (approximately)
elapsed = preferred_clock() - start

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/requests/adapters.py:696, in HTTPAdapter.send(self, request, stream, timeout, verify, cert, proxies)
   resolved_timeout = TimeoutSauce(connect=timeout, read=timeout)
try:
--> 696     resp = conn.urlopen(
       method=request.method,
       url=url,
       body=request.body,  # type: ignore[arg-type]  # urllib3 stubs don't accept Iterable[bytes | str]
       headers=request.headers,  # type: ignore[arg-type]  # urllib3#3072
       redirect=False,
       assert_same_host=False,
       preload_content=False,
       decode_content=False,
       retries=self.max_retries,
       timeout=resolved_timeout,
       chunked=chunked,
   )
except (ProtocolError, OSError) as err:
   raise ConnectionError(err, request=request)

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/urllib3/connectionpool.py:793, in HTTPConnectionPool.urlopen(self, method, url, body, headers, retries, redirect, assert_same_host, timeout, pool_timeout, release_conn, chunked, body_pos, preload_content, decode_content, **response_kw)
response_conn = conn if not release_conn else None
# Make the request on the HTTPConnection object
--> 793 response = self._make_request(
   conn,
   method,
   url,
   timeout=timeout_obj,
   body=body,
   headers=headers,
   chunked=chunked,
   retries=retries,
   response_conn=response_conn,
   preload_content=preload_content,
   decode_content=decode_content,
   **response_kw,
)
# Everything went great!
clean_exit = True

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/urllib3/connectionpool.py:470, in HTTPConnectionPool._make_request(self, conn, method, url, body, headers, retries, timeout, chunked, response_conn, preload_content, decode_content, enforce_content_length)
try:
   # Trigger any extra validation we need to do.
   try:
--> 470         self._validate_conn(conn)
   except (SocketTimeout, BaseSSLError) as e:
       self._raise_timeout(err=e, url=url, timeout_value=conn.timeout)

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/urllib3/connectionpool.py:1125, in HTTPSConnectionPool._validate_conn(self, conn)
# Force connect early to allow us to validate the connection.
if conn.is_closed:
-> 1125     conn.connect()
# TODO revise this, see https://github.com/urllib3/urllib3/issues/2791
if not conn.is_verified and not conn.proxy_is_verified:

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/urllib3/connection.py:827, in HTTPSConnection.connect(self)
try:
   sock: socket.socket | ssl.SSLSocket
--> 827     self.sock = sock = self._new_conn()
   server_hostname: str = self.host
   tls_in_tls = False

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/urllib3/connection.py:239, in HTTPConnection._new_conn(self)
"""Establish a socket connection and set nodelay settings on it.

:return: New socket connection.
"""
try:
--> 239     sock = connection.create_connection(
       (self._dns_host, self.port),
       self.timeout,
       source_address=self.source_address,
       socket_options=self.socket_options,
   )
except socket.gaierror as e:
   raise NameResolutionError(self.host, self, e) from e

File /opt/hostedtoolcache/Python/3.12.15/x64/lib/python3.12/site-packages/urllib3/util/connection.py:73, in create_connection(address, timeout, source_address, socket_options)
if source_address:
   sock.bind(source_address)
---> 73 sock.connect(sa)
# Break explicitly a reference cycle
err = None

KeyboardInterrupt: 

Display Metadata¶

Let’s examine the simulation parameters:

print("Simulation Metadata:")
print('=' * 60)
for k, v in icc.metadata.items():
    print(f'{k:20s} : {v}')
print('=' * 60)

Visualize the Waveform¶

Plot the (2,2) strain and Ψ₄ modes stored in the catalog:

# Merger time from the stored strain
tmrg, _, _, _ = icc.find_max(kind='global')

fig, axes = plt.subplots(2, 1, figsize=(12, 8), sharex=True)

axes[0].plot(icc.u - tmrg, icc.hlm[(2, 2)]['real'], label='Re', linewidth=1.5)
axes[0].plot(icc.u - tmrg, icc.hlm[(2, 2)]['A'], label='Amplitude', linewidth=1.5, c='k', ls='--')
axes[0].set_ylabel(r'$h_{22}$', fontsize=14)
axes[0].legend(fontsize=11)
axes[0].grid(True, alpha=0.3)

axes[1].plot(icc.u - tmrg, icc.psi4lm[(2, 2)]['real'], label='Re', linewidth=1.5, c='C1')
axes[1].plot(icc.u - tmrg, icc.psi4lm[(2, 2)]['A'], label='Amplitude', linewidth=1.5, c='k', ls='--')
axes[1].set_ylabel(r'$\psi_4^{22}$', fontsize=14)
axes[1].set_xlabel(r'$t-t_{\rm mrg}$ (M)', fontsize=14)
axes[1].legend(fontsize=11)
axes[1].grid(True, alpha=0.3)

fig.suptitle(f'ICCUB Simulation ID: {ID:04d}', fontsize=18)
plt.tight_layout()
plt.show()

Testing the Integration Method¶

The catalog also provides the raw Ψ₄ multipoles, so we can run PyART’s IntegrateMultipole (fixed-frequency integration) ourselves and inspect its outputs:

# Get Psi4 data
t = icc.t_psi4
psi4 = icc.psi4lm[(2, 2)]['z']
l, m = 2, 2

integr_opts = {
    'f0': 0.002,
    'extrap_psi4': False,  # the catalog data is already extrapolated to r=infinity
    'method': 'FFI',
    'window': [20, -20],
}

# Run the integration explicitly
mode = IntegrateMultipole(
    l, m, t, psi4,
    **integr_opts,
    mass=icc.metadata['M'],
    radius=1.0,
    integrand='psi4'
)

# Plot the integrated quantities
plt.figure(figsize=(14, 10))

plt.subplot(3, 1, 1)
plt.plot(t, mode.h.real, c='b', linewidth=2, label='Integrated')
plt.plot(t, icc.hlm[(2, 2)]['real'], label='Original', linewidth=1.5, color='r', linestyle='--')
plt.ylabel(r'Re[$h_{22}$]', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Strain', fontsize=14)
plt.legend()

plt.subplot(3, 1, 2)
plt.plot(t, mode.dh.real, c='b', linewidth=2)
plt.ylabel(r'Re[$\dot{h}_{22}$]', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Strain Derivative', fontsize=14)

plt.subplot(3, 1, 3)
plt.plot(t, mode.psi4.real, c='b', linewidth=2, label='Integrated')
plt.plot(t, icc.psi4lm[(2, 2)]['real'], c='r', linewidth=1.5, linestyle='--', label='Original')
plt.xlabel('Time (M)', fontsize=14)
plt.ylabel(r'Re[$\psi_4^{22}$]', fontsize=14)
plt.grid(True, alpha=0.3)
plt.title('Weyl Scalar', fontsize=14)
plt.legend()

plt.tight_layout()
plt.show()

The choice of \(f_0\) in the integration options determines the low-frequency cutoff for the waveform integration. A smaller value of \(f_0\) allows for more accurate integration of the low-frequency components, but may also introduce numerical noise. For instance, the real part of the integrated strain does not exactly taper to zero after merger. Choosing a slightly larger value of \(f_0\) can help mitigate this issue by suppressing the low-frequency noise

# Get Psi4 data
t = icc.t_psi4
psi4 = icc.psi4lm[(2, 2)]['z']
l, m = 2, 2

# Plot the integrated quantities
fig, ax = plt.subplots(figsize=(14, 6))
for i, f0 in enumerate([0.001, 0.002, 0.005, 0.01]):
    integr_opts['f0'] = f0

    # Run the integration explicitly
    mode = IntegrateMultipole(
        l, m, t, psi4,
        **integr_opts,
        mass=icc.metadata['M'],
        radius=1.0,
        integrand='psi4'
    )

    ax.plot(t, mode.h.real, linewidth=2, label=rf'$f_0={f0}$')
ax.set_xlabel('Time (M)')
ax.set_ylabel(r'Re[$h_{22}$]')
ax.grid(True, alpha=0.3)
ax.legend()

plt.show()

Summary¶

This tutorial demonstrated:

Downloading waveforms from the public ICCUB catalog
Integrating Ψ₄ (Weyl scalar) to obtain gravitational wave strain
Comparing the stored strain with the raw Ψ₄ data
Running PyART’s integration tools explicitly

Key Points¶

The public ICCUB release currently provides multipoles up to l=4
The Fixed Frequency Integration (FFI) method provides accurate strain from Ψ₄
Integration requires careful choices of f₀ (low-frequency cutoff) and time windows
PyART’s integration tools handle extrapolation and windowing automatically

Next Steps¶

Explore different integration methods and parameters
Compare ICCUB waveforms with other catalogs
Perform mismatch calculations with EOB models