WHFast 2-body simulation runs slower than its N-body equivalent

[rcdorsey]

Environment

• REBOUND Version: 4.4.6
• API interface: Python
• Operating System: Linux, Ubunti 24.04.3 LTS

Describe the bug
A simulation with the Sun and test particles takes longer to integrate with WHFast than an N-body simulation (Sun, 8 planets and test particles). Varying the integration length, output timestep size, integration timestep size, and turning safe_mode off did not make the 2-body simulation faster than the N-body simulation.

To Reproduce
Python script: test_NEO_evolution_2vN.py
Data file: Granvik+_2018_Icarus.txt

Name/Affiliation
Rosemary Dorsey, University of Helsinki

[hannorein]

I'm not sure if the script is running as intended: it only takes 0.4s?

[rcdorsey]

Yes, with the current simulation parameters the simulation only takes 0.38s for a "2-body" simulation but the issue is that the analogous N-body simulation takes 0.36s. I am currently in a testing phase, and the aim is to integrate a large number of test particles (10^6 at a minimum) for a large number of timesteps (~10^6), so the runtime difference will become considerable, and I do not understand why the N-body simulation integrates faster.

[hannorein]

Thanks. It's usually better to benchmark simulations that run for longer. You don't want to be dominated by the setup and input/output. But I'll have another look...

[rcdorsey]

I have tried a variety of parameters, but all show the N-body simulation as faster

[hannorein]

Ok there are a couple of issue with your code.

1. Units: Right now you are units where G=1. Your code units for length and time are 1 AU and 1/(2pi) years respectively. Yet you choose a timestep of 10 which corresponds to about 1.59 years which corresponds to 7 steps per Jupiter orbit. That is not small enough. You need to have at least ~20 steps per smallest orbital period (which are your NEO, not Jupiter). I assume you mixed up days with year/2pi. Check out https://rebound.hanno-rein.de/units/
2. For the two-body integration, you should set your timestep equal to t_wanted - t_now so that you reach the next output time t_wanted with one step.
3. Don't use safe_mode=False unless you know what you're doing. Right now your outputs are not correct because you do not synchronize the simulation beforehand. The easiest way is to keep the safe_mode enabled.
4. When I run your test, at least 50% of the time is overhead, e.g. input, output, calculating orbital elements, etc. When benchmarking, always try to find a test where you spend most of the time in the routine you want to optimize.

I hope this helps.