Initial orbital parameters and eccentricity control

[nilsvu]

Currently we use some old SpEC scripts to do basic eccentricity control, both for setting the initial parameters and for the iterations. Specifically, we use ZeroEccParamsFromPN.py for the initial parameters (#5933) and old eccentricity fits for the iterations (#5895). We have several options for improving these:

• Use the eccentricity control based on EOB models that has been tested in SpEC.
• Use the new varpro fits that have been tested in SpEC.
• Extend the PN formulae to more terms and to include eccentricity.
• Extend either of these methods to work better for BNS/BHNS as well.

We probably want a combination of these methods, until it's clear that one works much better than the others.

These improvements are orthogonal to the eccentricity control automation procedure, which is tracked in #5909.

[nilsvu]

Ping @PunJustice @haraldp271 @vtommasini @shabibti @taknapp @kchatziioannou @guilara since you've all been involved in eccentricity control. Let me know if you have any input :)

[knelli2]

@markscheel as well just so you're aware

[nilsvu]

In a meeting today with @markscheel @nilsdeppe @knelli2 we came up with an interface that can be used for eccentricity control. We think that we can use any eccentricity control method as long as it conforms to this interface. We can also import Python code that's merged into SpEC as long as it conforms to this interface. Just note that SpEC is private, so code merged into SpEC can only be used within the collaboration, whereas code merged into SpECTRE is publicly available.

Just two functions are needed:

class OrbitalParameters:
    """These parameters define a binary orbit.

    They can be used to generate initial data and evolve it.
    Object A is initially at positive x, and object B is at negative x.

    Of the orbital parameters, only a subset may be specified. It is the responsibility
    of the functions below to check for consistency of these parameters
    and fill in the rest.
    """
    # Intrinsic parameters
    mass_ratio: float  # q = M_A / M_B >= 1
    dimensionless_spin_a: Sequence[float]
    dimensionless_spin_b: Sequence[float]
    # Orbital parameters
    separation: Optional[float]  # D_0
    orbital_angular_velocity: Optional[float]  # Omega_0
    radial_expansion_velocity: Optional[float]  # adot_0
    eccentricity: Optional[float]
    mean_anomaly_fraction: Optional[float]  # The mean anomaly divided by 2 pi, so a value between 0 and 1.
    num_orbits: Optional[float]
    time_to_merger: Optional[float]

initial_orbital_parameters(params: dict) -> dict
    """Set the initial orbital parameters for a binary simulation.

    Arguments:
      params: Dictionary of requested orbital parameters. Fill in the rest and return the dict.
    """

update_orbital_parameters(params: dict, trajectories: np.ndarray) -> dict
    """Update the orbital parameters based on simulation output.

    Arguments:
      params: Dictionary of current orbital parameters. Update them and return the dict.
      trajectories: 2D array with 7 columns: time, xyz coordinates of object A, xyz coordinates of object B.
    """

Note: this interface is probably going to change a bit once we get more experience. For example, we might include parameters for neutron stars and more simulation output data.

Some important requirements on these functions:

1. They must not print anything unless through logging.
2. They must not plot anything or write any files unless given options to do so.
3. They may depend on other Python packages, but consider the dependencies carefully and document clearly what they are and why you need them. We can fairly easily pull in pure-Python dependencies through pip, but may have trouble installing/compiling non-Python dependencies on supercomputers.
4. For code in SpEC, try to avoid importing lots of SpEC-dependent code. Rather, have a Python file that only provides this interface, and use it with SpEC in a different file.

Note that we're not asking anyone to implement this for us right now, just note that we can use any eccentricity control method in SpECTRE if this interface is provided. And we'd be happy to! So if you want to make your eccentricity control method available to SpECTRE, this interface is how to do it. You can merge it into SpEC, SpECTRE, or both.

[PunJustice]

This all looks great! One note is that if you are thinking of using waveform-based eccentricity control (the EOB based method I've recently implemented in SpEC), we will need access to the waveform instead of the trajectories. In SpEC we use the finite radius strain at the closest extraction sphere, however we plan on testing different extraction spheres/running CCE on partial waveforms in the future. All of this is to say that we wouldn't be able to interface with update_orbital_parameters as is, and we would need more simulation output data.

[nilsvu]

Ok I think we have to dispatch to the different methods anyway, so we can pass different arguments to them. I've also been thinking about the method where we run a grid of simulations in parallel to determine the next update from a finite difference approximation in parameter space, which needs input from multiple simulations.

So all the requirements listed above still hold, and the interface should be close to the one outlined above, but each method may have to deviate slightly in the arguments it takes.

[guilara]

I think this looks good. Thank you for the update. I need think a bit more, but for the method I was working with, I'm think I might need more params to manually specify the time window for the fits (but anyway, the ultimate goal would be to have this automated).