hannorein / rebound

💫 An open-source multi-purpose N-body code.
https://rebound.readthedocs.io/
GNU General Public License v3.0
849 stars 219 forks source link

Wrong binary load for z-axis #654

Closed malsuar closed 1 year ago

malsuar commented 1 year ago

Description: I initialised and simulated a system of particles and then saved the simulation in a binary file. I then loaded the binary into another program, but the initial conditions do not match in the z-direction. Both figures were generated by the same code.

These are the figures: rebound_original_report

rebound_load_report

hannorein commented 1 year ago

Hi Mario!

This is strange. Can you post some code along with this. Otherwise it's really hard for me to help you debug this.

malsuar commented 1 year ago

Of course, I can. Let me simplify the code a bit. Thanks!

malsuar commented 1 year ago

In the meantime, this is basically how I load the simulation:

t_max = tiempo
sa       = rebound.SimulationArchive("./bin/bkp3/chkpd_%s_%d.bin"%(moon.name,t_max))
sim     = sa.getSimulation(t=0)

for i in range(ndp):
    xprm[i] = (sim.particles[i+2].x-sim.particles[0].x)/moon.radius 
    yprm[i] = (sim.particles[i+2].y-sim.particles[0].y)/moon.radius 
    zprm[i] = (sim.particles[i+2].z-sim.particles[0].z)/moon.radius 

rprm = (xprm**2+yprm**2+zprm**2)**0.5   

And finally, I create the figure in this way:

plt.figure(figsize=(12,4),constrained_layout = True)

ax1=plt.subplot(131)
ax1.scatter(rprm,xprm, s=0.4)
ax1.set_ylabel("x")
ax1.set_xlabel("r")

ax2=plt.subplot(132)
plt.title("Loaded REBOUND at t=0")
ax2.scatter(rprm,yprm, s=0.4)
ax2.set_ylabel("y")
ax2.set_xlabel("r")

ax3=plt.subplot(133)
ax3.scatter(rprm,zprm, s=0.4)

ax3.set_ylabel("z")
ax3.set_xlabel("r")

plt.savefig('./figures/rebound_load_report.png')
malsuar commented 1 year ago

The content of zprm 'is' zero:

image

hannorein commented 1 year ago

That looks all good. So the issue must be somewhere else.

hannorein commented 1 year ago

I'm pretty sure this is not a bug related to binary files. The particles probably had already a z position of 0 before you've saved the simulation to a file.

malsuar commented 1 year ago

How odd!

For instance, this is how I save the binary in the original program:

sim.automateSimulationArchive("./bin/chkpd_%s_%d.bin"%(moon.name,t_max),interval=0.1,deletefile=True)

And this is the plot created at the end of the simulation, for t=100. Note that z !=0: Final_configuration_301_b

On the other hand, this is what I get, after reading the binary at t=100. detail_Moon_t60

Also, snapshots for t > 0 have values for z != 0. for t=50 i get:

detail_Moon_t50

hannorein commented 1 year ago

I'm not sure what could be the issue. If I just save something to a file and read it back in, it works:

import rebound
sim = rebound.Simulation()
sim.add(x=1,y=2,z=3)
sim.save("test.bin")
sim2 = rebound.Simulation("test.bin")
sim2.status()

similarly:

import rebound
sim = rebound.Simulation()
sim.add(m=1)
sim.add(a=1,inc=0.1,omega=0.1)
sim.automateSimulationArchive("test.bin",interval=1.,deletefile=True)
sim.integrate(10)

sa = rebound.SimulationArchive("test.bin")
sim2 = sa.getSimulation(t=0)
sim2.status()

So I'm afraid I need enough of your code to reproduce the issue myself.

malsuar commented 1 year ago

Mmmm... These lines also work well for me. Even when I change the integrator, or get the status from the Horizon system. I will create a more refined version of the code to share with you. Thank you so much!

hannorein commented 1 year ago

Any update?

malsuar commented 1 year ago

Hi, Hanno. You were right from the beginning! I tried loading the binary files from other simulations of the same dataset and it worked quite well.
I also ran the simulation again, and it worked without problems. However, I couldn't figure out what was wrong with that first experiment.

Thank you very much for your time and cooperation. I will let you know if anything strange happens. ':)

hannorein commented 1 year ago

Glad to hear it was resolved!