What is the best way to migrate part of a Simulationarchive?

[astroboylrx]

Environment
Which version of REBOUND are you using and on what operating system?

• REBOUND Version: 4.4.3
• API interface: Python
• Operating System (including version): Linux (WSL+Ubuntu 22.04, HPC+Rocky Linux release 8.10), MacOS 15.1

---

I have a Simulationarchive that I want to extract the first N snapshots and save them to a new Simulationarchive (b/c I messed up the remaining snapshots). What would be the best way to do so?

The original Simulationarchive was created by a few simulation runs, each with regular intervals. I tried the most naive way:

sa = rebound.Simulationarchive("data.bin")
for ids in range(N):
    sa[ids].save_to_file("test.bin")

Then I went to verify if I got them right

foo = rebound.Simulationarchive("test.bin")
for ids in range(0, N):
    if not sa[ids] == foo[ids]:
        print(ids)

Weirdly, the first M snapshots are the same, and then the rest of them are different. The exact value of M seems to be somewhat random (I haven't figured out the pattern). And here is what diff looks like (similar for all snapshots later than M)

>>> sa[980].diff(foo[980])
ri_ias15.at:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.x0:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.v0:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.a0:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.csx:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.csv:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.csa0:
< (168 bytes, values not printed)
---
> (168 bytes, values not printed)
ri_ias15.g:
< (1176 bytes, values not printed)
---
> (1176 bytes, values not printed)
ri_ias15.b:
< (1176 bytes, values not printed)
---
> (1176 bytes, values not printed)
ri_ias15.csb:
< (1176 bytes, values not printed)
---
> (1176 bytes, values not printed)
ri_ias15.e:
< (1176 bytes, values not printed)
---
> (1176 bytes, values not printed)
ri_ias15.br:
< (1176 bytes, values not printed)
---
> (1176 bytes, values not printed)
ri_ias15.er:
< (1176 bytes, values not printed)
---
> (1176 bytes, values not printed)

They looks like something related to IAS15 but I'm using mercurius.

I then tried

sa[999].save_to_file("./test2.bin")
dan = rebound.Simulationarchive("./test2.bin")

It turned out sa[999] is still different than dan[0], and sa[999].diff(dan[0]) printed very similar output to above.
Only when I tried

dan[0].save_to_file("./test3.bin")
egg = rebound.Simulationarchive("./test3.bin")

I got dan[0] == egg[0] being True.

I also noticed that, with the same number of snapshots, test.bin has a larger file size than data.bin. Why is that?

Thank you very much in advance!

[hannorein]

Hm. What you're doing should work. I'm not quite sure what's going on. I'll need to investigate...

(FWIW Mercurius uses IAS15 internally.)

[astroboylrx]

I see. Thanks for the explanation! (Sorry I should read the doc better)

In case it helps, I put two example Simulationarchive files here

• in the first data file, the first different snapshot happens at index 980;
• in the second file, first different snapshot happens at index 4977.

[hannorein]

Ok. I can reproduce this:

import rebound
sim = rebound.Simulation()
sim.add("outer solar system")
for i in range(1, sim.N):
    sim.particles[i].m = 2e-2 # make sure IAS15 is triggered
sim.integrator = "mercurius"
sim.dt = sim.particles[1].P/20
sim.save_to_file("test.bin", step=1, delete_file=True)
sim.integrate(100)

sa = rebound.Simulationarchive("test.bin")
for ids in range(len(sa)):
    sa[ids].save_to_file("test2.bin")
foo = rebound.Simulationarchive("test2.bin")
for ids in range(0, len(sa)):
    if not sa[ids] == foo[ids]:
        print(ids)

[hannorein]

Still no solution, but I've narrowed it down a bit. Strangely, it only occurs when using mercurius, not when using IAS15 by itself. The results are often not reproducible, hence this is probably a memory corruption issue. Also, the same issue occurs when running something as simple as (no copying of any simulationarchive required):

sa = rebound.Simulationarchive("test.bin")
for ids in range(0, len(sa)):
    a = sa[ids]
    b = sa[ids]
    if a != b:
        print(ids)

[hannorein]

Shortest example so far:

sim = rebound.Simulation()
sim.add("outer solar system")
for i in range(1, sim.N):
    sim.particles[i].m = 2e-2
sim.integrator = "mercurius"
sim.step()
sim.save_to_file("test.bin", delete_file=True)
sa = rebound.Simulationarchive("test.bin")
a = sa[0]
b = sa[0]
print(a.diff(b))

[astroboylrx]

For this shortest example, if I run one step further, then they seem to converge.

a = sa[0]
a.step()
b = sa[0]
b.step()
a.diff(b)

output:

walltime:
< 1.738000e-03
---
> 1.837000e-03

Similar things happen for my previous examples, though much more steps are needed (one needs to run 1000 years and the other needs 100 years).
If the differences can disappear, does that mean those ri_ias15.XXX won't affect orbital integration results?

Regarding "memory corruption issue", does that have anything to do with the different file sizes? Simply iterating through a Simulationarchive and save_to_file yield a new Simulationarchive that is about twice larger in size.

[hannorein]

The IAS15 routines are reset at the beginning of the MERCURIUS step, so the differences might not affect the simulation when you continue. But if the memory is corrupted, then anything can happen...

[hannorein]

It took me a while to find the bug but it should be fixed now. Thanks for reporting it!

[astroboylrx]

Many thanks for fixing it!