Include "POWER" code to extrapolate waveforms to scri+ into ET

[Roland Haas]

POWER is an open source, python package to monitor the status and progress of numerical relativity simulations, and to post-process the data products of these simulations to compute the gravitational wave strain at future null infinity.

It is available from its git repository at NCSA: https://git.ncsa.illinois.edu/elihu/Gravitational_Waveform_Extractor/-/tree/rhaas/POWER2

There are two ongoing enhancements for it undeway that may be included in a final proposal if they finish quickly enough:

• full conversion to LIGO NR Waveform format: https://arxiv.org/abs/1703.01076
• offer method by Nakano et al. https://arxiv.org/abs/1503.00718 as an alternative to power-law extrapolation to scrI+.

‌

[Roland Haas]

Zach agreed to review POWER.

[Roland Haas]

Hold until there is a cleaned up version.

[Roland Haas]

The code is currently in a branch rhaas/POWER2 which you can check out like so:

git clone -b rhaas/POWER2 https://git.ncsa.illinois.edu/elihu/Gravitational_Waveform_Extractor

To test it you can use the waveform data stored in the Zenodo repo for the ET GW150914 example:

curl 'https://zenodo.org/record/155394/files/GW150914_28.tar.bz2?download=1' | \
  tar -vjx GW150914_28/output-000{0..4}/GW150914_28/{mp_psi4.h5,quasilocalmeasures-qlm_scalars..asc}  GW150914_28/output-0000/GW150914_28/TwoPunctures.bbh GW150914_28/SIMFACTORY/properties.ini

which extracts just the waveform, QuasiLocalMeasures and TwoPunctures files.

Once the data is there your can run the power-law extrapolation using:

./power.py --radii '[2:5]' --method POWER GW150914_28

which will produce and output directory Extrapolated_Strain. This can be compared to results obtained eg from SimulationTools: https://simulationtools.org/documentation/english/tutorials/waveforms :

<<SimulationTools`
$SimulationPath = {"."};

SimulationOverview["GW150914_28"]

extrapolationorder = 2;
omega0 = 0.025;

hExt = ReadRadiallyExtrapolatedStrain["GW150914_28", 2,2, omega0, extrapolationorder, Radii -> {136.0, 167.0, 214.0, 300.}]

ListLinePlot[Re[{hExt}]]

The sample output wold look like this from gnuplot + POWER

and from SimulationTools:

‌

[Roland Haas]

Missing things:

1. full documentation
2. regression test
3. documented test against SimulationTools
4. typos and other things in --help output
5. check list of authors, paper describing original POWER version

‌

[Maria]

Hi Roland,

Thank you for sending me the link. I read the paper arXiv:1708.0294, and download the Power using those instructions too. That version of Power has eccentric data, and the paper also makes reference to the data for a test BNS merger with a hybrid equation of state (MS1b).

I've been working for the past 2 weeks to add several hybrid EOSs to the parameter file provided of the Gallery web site, by using the parameter files that come with the Parma Code, and IllinoisGRMHD. It would be great if I would have this par file for comparison as well. Will you please give me access to that parameter file?

thanks, Maria

PS: the Power from the paper also has plot.py and monitor_waveform.sh. I assume those are not completely ready to be included in the new ETK release? I assumed that plot.py is important. I get an error with monitor_waveform.sh.

./monitor_waveform.sh simulations/J0040_N40

./monitor_waveform.sh: line 3: 25330 Segmentation fault: 11 ./plot_monitor.py $run

[Roland Haas]

I removed those two since they are not quite related to POWER and were very specific to the NCSA’s simulations, though getting a SEGFAULT is not something I would have expected.

monitor_waveform.sh in particular did not work as intended since the fixed-frequency-integration used in POWER does not work well for simulations that have not yet finished (since it uses no windowing so if a signal does not start and end with 0 amplitude then the FFT produces artifacts).

[Roland Haas]

Would you think that code is good enough that we can include it in the ET”s master for people to test (note that this does not constitute any statement about the code be ready to be included in the ET release, just that including it in master “does no harm”).

[Maria]

A comparison with Simulation Tools would be useful. I’ll be working on that. As for documentation, what you included in this ticket could be used. Is is short, but has the essentials.

I do have the following questions to begin with:

1. when pulling the data for GW150914_28 there are 5 output-000* My question is: how do they differ? Is psi4 extracted at a different radius in each of those runs? Could this be made clearer?
2. the data pulled included, besides mp_psi4.h5, TwoPunctures.bbh and quasilocalmeasures-qlm_scalars..asc. Are those necessary?
3. the path to ml_psi4.h5 is hard coded as GW150914_28/output-000*/GW150914_28/. Most likely users will want flexibility with the name of the directories in which ml_psi4.h5 is located. Can this be changed? Or, if not, instruct the user to bring the ml_psi4.h5 used for extrapolation in a specific directory?
4. Could a snippet of the output that needs to be included in the parameter file to generate mp_psi4.h5 be included?
5. Would it be possible to extrapolate mp_psi4 in ascii format?

I do hope those questions are easy to answer.

[Roland Haas]

Thanks for the comments. Some of them are easy to address, some are not quite so straightforward.

1. each output-???? directory corresponds to one segment of the simulation for this simulation which took 5 segments to complete. So each one contains a different time-range of the data with 0000 being the oldest and 0004 being the last and final part. These are the raw simulation output and the different extraction radii for psi4 are different datasets in the HDF5 file mp_psi4.h5. Definitely this should be documented as “Simfactory like layout” at the very least.
2. TwoPunctures.bbh is used by the code to obtain a measure of the total mass of the spacetime so that it properly rescale all data so that a total mass of 1 is used in the output files. Using the POWER method only requires TwoPunctures.bbh while Nakano’s method, since it needs information about the spin of the final black holes, requires quasilocalmeasures-scalars..asc as well.
3. the path to mp_psi4.h5 should not be hardcoded. The source code (see https://git.ncsa.illinois.edu/elihu/Gravitational_Waveform_Extractor/-/blob/rhaas/POWER2/power.py#L381) indicates that all files that are found by the shell glob output-????/*/mp_[pP]si4.h5 will be accepted. Older versions (in particular the ‘master’ branch) might have hardcoded this. If this still is hardcoded for you, can you point me to the source code line so that it can be fixed, please?
4. that is a very good idea indeed. Yes, this will be included in the documentation.
5. it would be possible but has not yet been implemented right now. Implementing that is probably not too hard since the files are only accessed in a handful of places ([1], [2], [3], [4], [5], [6]) that could easily be turned into functions that would accept both HDF5 or ASCII files.

‌

[Maria]

Roland,

You’re right, the path is not hard coded in power.py.

I was able to install Simulation Tools as well and check the Power extrapolated plot against it. You give

hExt = ReadRadiallyExtrapolatedStrain["GW150914_28", 2,2, omega0, extrapolationorder, Radii -> {136.0, 167.0, 214.0, 300.}]

I see there is a little shift in time. Can the command above be tweaked to improve this? Below is the plot.

Would be great if you can include the “plot.py” script as well in the Gravitational_Waveform_Extractor.

As for “monitor_waveform.sh”, if you prefer not to include it, would it be appropriate to cut “to monitor the status and progress of numerical relativity simulations“ from the description?

[Maria]

Roland, regarding your question:

“Would you think that code is good enough that we can include it in the ET”s master for people to test (note that this does not constitute any statement about the code be ready to be included in the ET release, just that including it in master “does no harm”).”

The answer is yes.

‌

[Roland Haas]

I had not even realized you had provided a plot (only looked at the email notification). Hmm, will have to see where the offset comes from. Physically it does not matter of course since it is just a time offset but the code was originally designed to reproduced SimulationTools results as much as possible.

[Maria]

Roland, I added the parameters for outputting psi4.h5 to the nsnstohmns.par. This is what I added:

______________________________
WeylScal4::fdOrder = 8
WeylScal4::calc_scalars = "psis"
WeylScal4::calc_invariants = "always"

Multipole::nradii = 7
Multipole::radius[0] = 100
Multipole::radius[1] = 115
Multipole::radius[2] = 136
Multipole::radius[3] = 167
Multipole::radius[4] = 214
Multipole::radius[5] = 300
Multipole::radius[6] = 500
Multipole::ntheta = 120
Multipole::nphi = 240
Multipole::variables = "WeylScal4::Psi4r{sw=-2 cmplx='WeylScal4::Psi4i' name='psi4'}"
Multipole::out_every = 512
Multipole::l_max = 4 #8
Multipole::output_hdf5 = yes

Multipole::output_ascii = no

___________________________

One other change I made was to increase the boundary points and ghost points to 4.

I do have a problem with the nsnstohmns.par grid settings. The way is it now: N=396, dx = 18

I am getting interpolation error:

WARNING level 1 from host babiuct7910-1 process 0
while executing schedule bin CCTK_ANALYSIS, routine Multipole::Multipole_Calc
in thorn AEILocalInterp, file /home/babiuc/Cactus/arrangements/Numerical/AEILocalInterp/src/Hermite/../template.c:1109:
->
CCTK_InterpLocalUniform():
interpolation point is either outside the grid,
or inside but too close to the grid boundary!
(this may be caused by a global interpolation with
driver::ghost_size too small)
0-origin interpolation point number pt=13469 of N_interp_points=13470
interpolation point (x,y,z)=(237.748,116.683,424.1)
grid x_min(delta_x)x_max = -36(9)423
grid y_min(delta_y)y_max = -423(9)423
grid z_min(delta_z)z_max = -36(9)423

WARNING level 1 from host babiuct7910-1 process 0
while executing schedule bin CCTK_ANALYSIS, routine Multipole::Multipole_Calc
in thorn Multipole, file /home/babiuc/Cactus/configs/bns/build/Multipole/interpolate.cc:10:
-> CCTK_InterpGridArrays returned error code -1002.

______________

However, if I change it to N=576, dx=18 the interpolation errors disappear. I got this number because I tried to run with a grid multiple of 2^n, for example N=512 and dx =16. However, when I do that, the stars origin is displaced from 15.1875 to 15.25. I do not understand why is that, but it seems important to run with a dx multiple of 6.

So, if we want to add Power to nsnstohmns.par, I suggest replacing N=396 to N=576.

One other improvement I suggest to nsnstohmns.par is to add horizon finder. I am going to run the par file and use Power on the generated psi4.h4.

Maria

‌

[Roland Haas]

I added a pull request for the thornlist change to:

https://bitbucket.org/einsteintoolkit/manifest/pull-requests/7/einsteintoolkitth-add-power-clean-up

which puts power.py into the utils directory and its “Analysis” subdirectory.

[Gabriele Bozzola]

I would like to use POWER to validate kuibit.

I downloaded the code and used it (but haven’t compared against kuibit yet) and I have some suggestions for future improvements:

• It would be nice to be able to control explicitly the parameters via command-line, e.g., the ADMmass or the cutoff frequency. Depending on TwoPunctures limits the applicability of the tool.
• The function getFinalSpinFromQLM assumes that scalars are output as one group per file
• If one passes as folder ., then the output will be hidden files
• It would be nice to be able use the code on arbitrary directory structures
• The code seems to make some assumptions (e.g, I see, junk_time = 50). It would be great if these were all clearly summarized somewhere. (For example, this junk time would lead to incorrect results if ADMmass = 0.1, but no error would be thrown.)

Other than that, thanks for the tool, I am looking forward to using it more carefully.

[Roland Haas]

I will take a look Sunday (busy with other stuff today).

[Roland Haas]

updated to contain:

• mp_psi4 in ascii format
• made TwoPucntures.bbh and QLM output optional for all methods by allowing user to specify ADM mass, cutoff frequency, final spin and mass explicitly in the Python interface
• added Python module option to set the glob used to find output directories, it defaults to output-????/*
• add Python module option and function arguments to set basename of psi4 files use, defaults to mp_[Pp]si4
• added code to handle simulation directories that are . and FOO/ ie end in a slash
• maked the junk time a module option

No docs yet :-(

‌

[Maria]

I tried to run power with nsnstohmns.

It is able to read *asc files, even defaults on them!

Unfortunately, it still complains about missing TwoPunctures.

‌

▶ ./power.py --radii '[2:5]' --method POWER nsnstohmns_dx9
Extrapolating with POWER method...
Traceback (most recent call last):
File "./power.py", line 860, in <module>
main()
File "./power.py", line 835, in main
strains = POWER(args.path, PSI4_GLOB, radii, modes)
File "./power.py", line 481, in POWER
meta_name = glob.glob(os.path.join(simdirs, "TwoPunctures.bbh"))[0]
IndexError: list index out of range

‌

Regarding the path, it still requests a nested path for the files of type: ‘*/output-????/*/’

It would be good to be able to receive just the mp_Psi4.h5.

The error is either:

It gives me either:
FileNotFoundError: [Errno 2] No such file or directory

or

error: argument path: invalid dir_path value

‌

Would be good to explain plain in the documentation about this.

[Gabriele Bozzola]

I also tried the new POWER. It took me a while to figure out how to set f0 and M_ADM. I would have expected to find these options among the ones you can specify via command-line, or at least among the “module options”. I saw that there are checks like f0 == FROM_TWOPUNCTURES, so I even considered whether I had to set f0 by setting the variable FROM_TWOPUNCTURES.

[Roland Haas]

I will have to write proper documentation still., which will document these,

I am be quite reluctant to add new options to the command line since it can easily become a mess. For more complex use cases, using it as a Python module and calling the POWER or Nakano routines directly is the preferred way of using it.

The thinking is that “module options” will be things that would be identical between simulations while function arguments are things that change from one simulation to the next.

In that respect setting f0 would be a command line option candidate (specific to fixed frequency integration, though changes from simulation to simulation).

[Gabriele Bozzola]

What I had in mind is that one can put power.py in their PATH and use it for all their simulations without having to touch the code. This requires the parameters to be adjustable via command-line. But if that is not the designed way of using the tool, that would be fine.

Also, for clarity, instead of defining a string FROM_TWOPUNCTURES that is going to be compared to something that is morally a float (f0 and ADMMass), I would recommend setting the default of two parameters to None in the function definitions, and then checking if value is None. If not, one can assume that the values from TwoPunctures are used (which is what is happening now). So, something along the lines of

if f0 is None:
  f0 = getCutoffFrequencyFromTwoPuncturesBBH(meta_name)

(As someone that is not familiar with the code, I found the FROM_TWOPUNCTURES variable very confusing, but maybe its role will be clarified by the documentation.)

Two additional comments are:

• When the spin is computed from QLM a specific column is read. In my simulation, column 59 corresponds to 59:qlm_coordspinx[1], which is not the mass. Some assumptions are being made here and no checks are performed. This might lead to subtly wrong results.
• It would be nice, and probably easy, to extend POWER to be able to compute the strain at a fixed radius (essentially, only doing the FFI). My use case is this. I tried POWER on my data and the output doesn’t look good at all. I don’t know if the problem is the interpolation to infinity or the data. Maybe the data becomes unreliable at a specific distance. If I can plot the strains at all the various distances I can get a sense of what is going on.

‌

[Roland Haas]

Maria, Miguel since you are the first non-local users, which options would you find useful on the command line? Clearly setting f0 seems among that list. Any others?

[Roland Haas]

Right now using FROM_TWOPUNCTURES serves the purpose that is being suggested for f0 = None (FROM_TWOPUNCTURES is just a string and I only use the variable name to try and avoid typos as much as possible given that Python is fully dynamically typed and requires no variable declarations at all). I would rather not use None for this since there is only one None and one can get the total mass also in other ways which is easy to do by having another FROM_XXX but not if None is already taken.

You can use power.py as a Python module import power then call the Python functions directly. This means you have to write Python code rather than shell code.

I have code around somewhere that would parse the header and then picks the correct column. That should be easy to incorporate.

My first guesses for bad POWER output would be:

• make sure your psi4 starts and ends at 0 amplitude
• make you may have to change the list of detectors being used, using detectors too far out or in can give bad results
• TODO: probably there should be options to control the order of the polynomials that are being fit to the data

‌

[Gabriele Bozzola]

Right now using FROM_TWOPUNCTURES serves the purpose that is being suggested for f0 = None (FROM_TWOPUNCTURES is just a string and I only use the variable name to try and avoid typos as much as possible given that Python is fully dynamically typed and requires no variable declarations at all). I would rather not use None for this since there is only one None and one can get the total mass also in other ways which is easy to do by having another FROM_XXX but not if None is already taken.

I understand the reasoning now and it makes sense.

You can use power.py as a Python module import power then call the Python functions directly. This means you have to write Python code rather than shell code.

My personal preference would be to call power with something like power.py …. --f0 0.1 --ADM 1 --method POWER . as opposed to creating a Python script, figuring out the function signature, taking care of saving the data myself, et cetera. Also, extending the main() function does not affect the functionalities for when you import power as a module, so you can still write Python code in exactly the same way.

‌

My first guesses for bad POWER output would be:

• make sure your psi4 starts and ends at 0 amplitude
• make you may have to change the list of detectors being used, using detectors too far out or in can give bad results
• TODO: probably there should be options to control the order of the polynomials that are being fit to the data

I’ll try to play around and see.

‌

[Maria]

If f0 is what will enable using Power to calculate the strain from BNS runs, then that’s it. Knowing the command line that should be used for this would be helpful too. And also, explaining how to manipulate Power in an IPython or Jupyter Notebook environment, making TwoPunctures.bbh and QLM output optional would be a good addition to the Documentation.

‌

[Roland Haas]

Documentation and command line options for overriding f0, M_ADM, M_final and a_final were added in git hash a39f1e8 "power: fix typo in docs" of Gravitational_Waveform_Extractor

[Maria]

Roland,

‌

In your documentation, you write: “Please consult the output of power.py --help for details on the supported command line options.“

Running ./power.py --help, as you instructed, gives me a long list of options, and I am overwhelmed. I read that I need to “Use
'QuasiLocalMeasures' or ‘QLM' to deduce fromquasilocalmeasures-qlm_scalars..asc (the default)”. I don’t have quasilocalmeasures-qlm_scalars..asc as output. Is there another way to recover those parameters without redoing the runs?

An example on how to use the script with BNS data will be much appreciated. If not, a short instruction on how do I find the steering parameters will be useful. Or, if this is possible only within Python, a short example of what functions should be called to obtain the necessary parameters.

‌

Oh, one more thing: please include in the documentation the essential lines that must be included in the par file, in order to be able to use POWER.

There is this:

______________________________
WeylScal4::fdOrder = 8
WeylScal4::calc_scalars = "psis"
WeylScal4::calc_invariants = "always"

Multipole::nradii = 7
Multipole::radius[0] = 100
Multipole::radius[1] = 115
Multipole::radius[2] = 136
Multipole::radius[3] = 167
Multipole::radius[4] = 214
Multipole::radius[5] = 300
Multipole::radius[6] = 500
Multipole::ntheta = 120
Multipole::nphi = 240
Multipole::variables = "WeylScal4::Psi4r{sw=-2 cmplx='WeylScal4::Psi4i' name='psi4'}"
Multipole::out_every = 512
Multipole::l_max = 4 #8
Multipole::output_hdf5 = yes

Multipole::output_ascii = no

___________________________

‌

Besides, shall one add quasilocalmeasures-qlm_scalars..asc?

[Roland Haas]

QLM output is only required for the Nakano method, not for the basic extrapolation method.

I added more documentation in git hash b23d519 "power: add some example parfile fragment" of Gravitational_Waveform_Extractor

[Roland Haas]

Checkout currently fails on Ubuntu / WSL with an invalid certificate error:

fatal: unable to access 'https://git.ncsa.illinois.edu/elihu/Gravitational_Waveform_Extractor.git/': server certificate verification failed. CAfile: none CRLfile: none

though curl https://git.ncsa.illinois.edu/elihu/Gravitational_Waveform_Extractor.git works fine (even if it only shows an HTML page stating that I will be redirected).