include PN expressions for initial momenta for BBH simulation

Issue #2417 resolved
Roland Haas created an issue

Post-Newtonian theory results in some of the longest and most complex mathematical expressions ever derived by humanity.

These expressions form the core of most gravitational wave data analysis codes, but generally the expressions are written in a format that is either inaccessible by others (i.e., closed-source) or written directly in e.g., C code and thus incompatible with symbolic algebra packages like Wolfram Mathematica or the Python-based SymPy.

Once in a symbolic algebra package, these expressions could be manipulated, extended, and output as more optimized C codes (e.g., using SymPy/NRPy+, thus speeding up gravitational wave data analysis).

This repository aims to provide a trusted source for validated Post-Newtonian expressions useful for gravitational wave astronomy, using the open-source SymPy computer algebra system, so that expressions can be output into Wolfram Mathematica or highly optimized C codes using the SymPy-based NRPy+.

Available from:

Comments (18)

  1. Barry Wardell

    The Black Hole Perturbation Toolkit has a similar repository within the context of pN expansions in the large mass ratio regime. The Post-Newtonian Self-force package is a Mathematica interface to the repository, but there is a desire within the community to have something which is tool-agnostic. There have also been discussions about extending the scope to cover more generic post-Newtonian series. It might make sense to coordinate those efforts with this one.

  2. Peter Diener

    I have reviewed the code. The code consists of a set of beautifully commented jupyter notebooks where, for each equation implemented, the reference to the paper where the equation was derived and presented is listed. Also where typos in published equations were found, this is clearly indicated. I have run the notebooks on my own laptop and didn’t find any discrepancies on any of the results (though I couldn’t get exporting to LaTeX and PDF to work; most likely to old software). As I see no reason why this code can not be included in the toolkit, I hereby approve it.

  3. Roland Haas reporter

    Looking a the information in this ticket (only, no other information having been provided) I am at this point missing (ie have no idea if they are present or not):

    • installation instructions
    • information on authors, license, copyright needed for inclusion in the ET and adding to the ET author list / ET announcement email
    • information on its dependencies (were those also reviewed?)

      • how are issues / tickets handled? Issue tracker in own repo / ET tracker?
    • a repository that can be tagged with the ET release tag and to which at least one ET maintainer (Zach I presume) has write access

    • the notebook seems to be an example but not actually documentation. It does not provide information for example what chi1U is or how the variables nchi2z etc. are being used. These may be present in the repo but where is not shown in the ticket.
    • tests on OSX

  4. Zach Etienne

    Roland, Peter, and Zach discussed this ticket during the October 15, 2020 ET weekly telecon. It was decided that Zach will need to

    • Modify GetComponents to clone the full NRPy+ repo, and create a symbolic link to the NRPyPN subdirectory to arrangements/EinsteinInitialdata
    • Copy the trivial NRPy+ functionality from the NRPy+ root directory to NRPyPN, so that NRPyPN can be fully standalone. Modify NRPyPN slightly so it no longer looks for that functionality in the parent directory.
    • Add a README to the NRPyPN root directory, containing author list; installation & usage instructions; and software license (2-clause BSD)
    • Discuss with Roland how to run Jupyter notebooks in Jenkins, as the Jupyter notebooks containing the expressions all contain extensive tests, which cause the Jupyter notebooks to error out if the tests fail. Note that these notebooks are already run within the NRPy+ Travis builds:

  5. Roland Haas reporter

    For just a subdirectory the thornlist fragment becomes very simple eg:

    # Component list for NRPyPN <>
    !CRL_VERSION = 1.0
    !DEFINE ROOT = Cactus
    !DEFINE ARR  = $ROOT/arrangements
    # Note: an alternative would be to include it in the "utils" directory in the
    # main Cactus directory since it is not really a thorn
    !TARGET   = $ARR/EinsteinInitialData
    !TYPE     = git
    !URL      =

    where the trick is to keep NRPyPN in the !CHECKOUT line so that it is not treated as thorn when compiling. This is similar to eg doc and par of McLachlan.

  6. Roland Haas reporter

    I looked a bit into what needs to be present. At least for my OSX+Homebrew test system I needed:

    • brew install jupyter
    • pip3 install sympy
    • sys.path.append(“/User/jenkins/Library/Python/3.8/lib/python/site-packages“) otherwise it would not find the package in jupyter notebooks (but would find in python scripts, very odd)
    • no other Python package seems required

    Python3.8 produced 2 syntax warnings of type ‘SyntaxWarning: "is" with a literatl. Did you mean "=="?’ in mpmath, one of them line 892. Probably not much that can be done about this though.

    The example notebook explicitly provides digits for the Euler-gamma constant. However sympy provides this as a symbol already sympy.S.EuerlGamma which probably could (and should) be used instead in the PN expressions (and Float(sympy.S.EuerlGamma) does give the correct digits).

    None of the Python functions seem to have a Python docstring, making eg help(pt.f_p_t) not very useful. is really used as some sort of include file injection a number of symbols (among them chi1U etc.) into the main namespace via from NRPyPN_shortcuts import * which is somewhat confusing.

    The pt module seems to have some internal state that is manipulated via things like pt.f_p_t. It would seem cleaner if instead there was a pt object that one created instead.

    There seems to be a typo in a comment in the example notebooks which reads:

    # Compute p_t, the radial component of momentum
    import PN_p_t as pt
    pt.f_p_t(m1,m2, chi1U,chi2U, r)
    # Compute p_r, the radial component of momentum
    import PN_p_r as pr
    pr.f_p_r(m1,m2, n12U,n21U, chi1U,chi2U, S1U,S2U, p1U,p2U, r)

    since it states that both p_t and p_r compute the radial component of the momentum. The first one, based on its name, is likely the tangential component instead.

  7. Roland Haas reporter

    @Zach Etienne added a file including:

    • license information
    • required sympy version
    • installation instructions
    • author list

  8. Roland Haas reporter

    Note that there is no (to my best knowledge) documentation in the sense that there is no document explaining (other than the example) how to use the script, what its offered functionality and its limitations are.

    There are no Python docstrings.

    There is no doc/documentation.tex which means NRPyPN will not show up on

  9. Zach Etienne

    @Roland Haas I like this idea. We hope mybinder will be around forever, but presumably users will notify us if/when there's a problem with mybinder.

  10. Log in to comment