Scientific Computing for Physicists  (Winter 2020) (Old site; new site is at https://scinet.courses)

For this assignment, you are to mimic one of the methods used by the Laser Interferometer Gravitational Wave Observatory (LIGO) for detecting gravitational waves (GWs) by analyzing signals from the detector, which led to the 2017 Nobel prize in physics for the detection of GWs.

We will use data generated by real numerical relativity simulations kindly made publicly available by the Center for Computational Relativity and Gravitation, see Campanelli, Lousto, Nakano, and Zlochower, Phys. Rev. D 79:084010, 2009).

An oversimplified representation of LIGOs "Matched Filtered Algorithm" goes as follows:

1. Consider a bank of predicted GW waveforms generated by numerical relativity simulations of binary black hole mergers. Each wave form is given as a complex-valued time series.
2. Collect a set of measured signals from the detector (also complex-valued time series).
3. Perform a power spectrum analysis, correlating the GWs signals with the signals in the simulated waveform bank. This correlation analysis is performed in the Fourier domain to make the results more robust, as the signals are usually embedded in noise.
4. Select candidate events with a high correlation coefficients with one of the predicted waveforms.

Assignment

In this assignment, you will implement something similar. You only need to consider one predicted GW signal from simulations of the merger of two black holes in the netcdf file GWprediction.nc, and 32 mock detector signals in the netcdf files detection01.nc, ..., detection32.nc. This data can be found on the teach cluster in the directory /scinet/course/phy1610/gwdata.

The file structure is as follows. Each file contains a variable called f which stores the values of the wave signal at regularly spaces time points. The file also contains a variable with the actual times at which the samples were taken but this information is not really needed in this assignment. The wave signal is a set of double-precision complex numbers, so each point has a double-precision real and imaginary part. Because NetCDF does not have native complex number support, the f variable is a nt x 2 matrix, where the time dimension nt is stored in the file. This is in fact also how complex numbers are laid out in memory in C++, so it's in fact fairly easy to read in these files, e.g., like this:

#include <rarray>
#include <netcdf>
#include <complex>
using namespace netCDF;
...
 NcFile ncfile("detection01.nc", NcFile::read);
 rvector<std::complex<double>> f(ncfile.getDim("nt").getSize());
 ncfile.getVar("f").getVar(f.data());
...

You will have to compute the correlations of the power spectrum of the prediction with each mock detection. The power spectrum F of a signal f is related to the fourier transform of that signal: for each wave number k (in this case they are really frequency numbers), the power spectrum is the square norm of the fourier component with that wave number, i.e.,

.

Note that the power spectrum is real-valued.

The correlation of two power spectra F and G is given by

where denotes an inner product, which here is simply,

Note that the denominator in the definition of is such that it is close to 1 for highly correlated F and G.

The algorithm you have to implement is structured as follows:

1. Read the predicted GW signal from GWprediction.nc.
2. Read one of the GW signal from observations detection01.nc . . . detection32.nc.
3. Compute the FFT of the two complex quantities, using FFTW.
4. Compute the power spectrum of both signals.
5. Compute the correlation coefficient between the power spectra as defined in Eq. (1), using a ?dot BLAS call for the inner product from Eq. (3).
6. Output the correlation coefficient
7. Repeat steps 2-to-6 for each of the signals in the observation set.
8. Finally, determine the 5 most significant candidates (those with the 5 largest values of the correlation coefficient) from the observations set.

For point 3, you will need an fft implementation. Usage the fftw module on the teach cluster for that.

For point 5, you will need a blas implementation. You should install openblas from source into (a subdirectory of) your $HOME directory on the teach cluster for that.

About the detected signals

The following information is not necessary for solving the assignment, but just in case you are interested, notice that the GW signal is a complex quantity, named , composed of , where and are quadrupole modes of the gravitational wave, i.e., of the perturbation of the metric describing the space-time far away from the source of the event. This event could, for instance, be the merger of two black holes.

What to submit when

As in previous assignments, you should still use git version control and a makefile for your code.

For simplicity, use one single git repository for the assignment.

To your git repository, add a plain text file with the results of your analysis, i.e. which are the 5 most significant candidates (those with the 5 largest values of the correlation coefficient) from the observations set.

Note that we expect you to use the best practices that we've been insisting on, like commenting your code and several meaningful git commit messages. While there is not a lot of opportunity for modularity here, stucturing your code into functions is a good idea.

Submit the "git2zip"-ed repo to the Dropbox by Monday March 9th, 2019 at 11:55 PM. The usual penalty for late submissions applies.

Last Modified: Friday Jun 5, 2020 - 14:54. Revision: 2. Release Date: Monday Mar 2, 2020 - 23:00.