Poor energy conservation with sink binary orbit in co-rotating frame

[themikelau]

I have simply been trying to evolve two sink particles in a ciruclar binary orbit. The energy conservation is really bad when in evolving in the co-rotating frame, with the binary out-spiralling within ~10 orbits and running into the energy conservation error. This does not happen in an inertial frame. I have identified that this issue originated from the pull request in https://github.com/danieljprice/phantom/commit/c9cc1eb89956d4f89c6807166be6c7cba6d8a7ab.

I have also tried setting use_fourthorder=.false. everywhere as per @danieljprice suggestion. But that did not resolve the problem, implying even the "false" case contained a bug.

[Yrisch]

Hi, Is there only two sinks in the simulation (just gravitationnal interaction) ? Do you have a setup to share ? It's strange that it don't work with use_fourthorder=.false.. It is supposed to be exactly the same as the previous implementation in that case. I will investigate...

[themikelau]

Hi @Yrisch , thanks for the reply. A simple way to reproduce this is to use SETUP=binary. I have attached the setup file and .in file in the zip file. I also had to fix a small bug along the way to get setup_binary to correctly set the orbit in the co-rotating frame, so please use the newly created binary branch.

In the setup, M1=M2=a=1, corotate=T. I have specified omega_corotate = 1.41421356 and iexternalforce = 2. In the corotating frame, the point masses should not move at all as they are in a circular orbit. But I observe significant deviations from around t ≈ 12 ≈ 6 hr. binary_sink_setup.zip

[Yrisch]

Thanks for the files ! I tried your setup and it also stopped but just after the first dump in my case... I supposed that the issue is in the ptmass_vdependent_correction routine. It's the only thing that has not been properly tested yet... I'm on it.

[themikelau]

I had to set export I_WILL_NOT_PUBLISH_CRAP=yes to keep the simulation running past the energy non-conservation threshold. Thanks for checking!

[Yrisch]

Okey I found the bug... It was actually in the ptmass_vdependent_correction. But I think that there was a more ancient bug into the code before my update. In this routine we compute force corrections due to coriolis. This correction seems to be added twice ( modifying the force by adding the corrections inside the routine but also a second time outside the routine). This doubled correction didn't annoy the previous version of the code because it used a local variable to immediately update the velocity of the sink. in the new version the correction is directly added into the force array. It then may propagate some errors into the code but I'm not completely sure about the exact reasons. Nevermind, by deleting the double correction I have a agreement between the previous and the current version of the code.

[Yrisch]

see #564 for more details in the code...

[themikelau]

Thanks for fixing this!

[danieljprice]

Thanks Yann for the speedy response here. I would mainly add that we need a unit test for this!!