Impossible to read more than one snapshot with python

[Francyrad]

Dear users

I tried to read the snapshots of my "archive.bin" with python, but i can only generate the picture of the first snapshot. When trying to reading the second. The code goes in segmentation fault. This problem is both in linux and MacOS. I have no idea if this is fault of my code. Can you please check?

import numpy as np
import rebound
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
import warnings

def main():
    archive_file = "archive.bin"
    num_snapshots = get_number_of_snapshots(archive_file)
    print(f"Number of snapshots: {num_snapshots}")

    for i in range(num_snapshots):
        print(f"Processing snapshot {i}")
        sim = get_simulation_at_snapshot(archive_file, i)
        plot_simulation(sim, i)
        print(f"Snapshot {i} processed")

def get_number_of_snapshots(archive_file):
    sa = rebound.SimulationArchive(archive_file)
    return len(sa)

def get_simulation_at_snapshot(archive_file, snapshot):
    sa = rebound.SimulationArchive(archive_file)
    sim = sa[snapshot]
    return sim

def plot_simulation(sim, snapshot_number):
    fig, ax = plt.subplots(figsize=(10, 10))

    R_earth_to_meters = 6371000
    ax_lim = 5 * R_earth_to_meters  # 5 Earth radii in meters

    ax.set_xlim(-ax_lim, ax_lim)
    ax.set_ylim(-ax_lim, ax_lim)
    ax.set_xlabel(r"x Position ($R_\oplus$)")
    ax.set_ylabel(r"y Position ($R_\oplus$)")

    for i, particle in enumerate(sim.particles):
        radius = particle.r
        if i == 0:  # Assuming the first particle is the planet
            color = 'white'
        else:
            color = 'gray'

        circle = plt.Circle((particle.x, particle.y), radius, color=color)
        ax.add_artist(circle)

    ax.set_title('Origin of the Moon. (C) Francesco Radica - Francyrad', color='w', fontsize=22)
    ax.set_facecolor("black")
    fig.patch.set_facecolor('black')

    ax.yaxis.label.set_color('w')
    ax.xaxis.label.set_color('w')
    ax.tick_params(axis='x', colors='w',labelsize=20)
    ax.tick_params(axis='y', colors='w',labelsize=20)
    ax.yaxis.set_major_formatter(ticker.FuncFormatter(lambda x, pos: f'{x / 1e7:.0f}'))
    plt.gca().set_aspect('equal', adjustable='box')

    plt.savefig(f'snapshot_{snapshot_number}.png', dpi=300, facecolor='black', bbox_inches='tight')
    plt.close(fig)

def empty_function(sim_pointer):
    pass

if __name__ == "__main__":
    main()

I got this error and the following warnings:

python3.10 snaps_pic.py            
/opt/homebrew/lib/python3.10/site-packages/rebound/simulationarchive.py:99: RuntimeWarning: Binary file was saved with a different version of REBOUND. Binary format might have changed.
  warnings.warn(message, RuntimeWarning)
Number of snapshots: 20
Processing snapshot 0
/opt/homebrew/lib/python3.10/site-packages/rebound/simulationarchive.py:146: RuntimeWarning: Binary file was saved with a different version of REBOUND. Binary format might have changed.
  warnings.warn(message, RuntimeWarning)
/opt/homebrew/lib/python3.10/site-packages/rebound/simulationarchive.py:146: RuntimeWarning: You have to reset function pointers after creating a reb_simulation struct with a binary file.
  warnings.warn(message, RuntimeWarning)
/Users/francyrad/Documents/rebound-main/examples/moon/snaps_pic.py:28: MatplotlibDeprecationWarning: The resize_event function was deprecated in Matplotlib 3.6 and will be removed two minor releases later. Use callbacks.process('resize_event', ResizeEvent(...)) instead.
  fig, ax = plt.subplots(figsize=(10, 10))
Snapshot 0 processed
Processing snapshot 1
zsh: segmentation fault  python3.10 snaps_pic.py

Thank you in advance for your support

[hannorein]

Hi Francesco,

I don't see anything obviously wrong with your code. I'm sure you've noticed the warning message about the version mismatch? What version of REBOUND was used to generate the archive file and what version are you using to read the fie in? If that's not the source of the problem, then you'd need to send me your binary file so I can run the code myself.

Hanno

[Francyrad]

I did 2 tries, both with rebound 3.24.2 that (should) be the last version available in python

Regarding the program, the code is been written in C. I did a simulation with the version of 3 April 2023 (i suppose it is the 3.24.2). Right now i downloaded the 3.24.3. I tried to run my code and i've the same problem.

If you can please try to read my archive.bin. I also share with you my code, in way that you can check if there is something wrong, the makefile and the python script

Thank for your availability and your help

Francesco

help.zip

[hannorein]

Thanks. The problem is that you have at least two particles which have exactly the same position and as a result the tree code ends up in an infinite loop until you reach a recursion limit or run out of memory. A better warning message would be nice, but the underlying problem has probably nothing to do with the Simulation Archive.

I can see in you code that you have implemented your own collision resolve function. I'm not quite sure what you're trying to do with it. I suspect there is some confusion about its use. Most importantly, the function has the wrong return type. It should return an integer that determines if one of the two particles should get removed from the simulation. You need to fix that to avoid unexpected results. Look at the implementation of the hard sphere routine that comes with REBOUND.

In the meantime, I'll try to add some error message and avoid the segmentation fault in this case...

[Francyrad]

Thank you! Solving the collision issue solved the problem! I close this! I wish you a good Job!

Francesco

[hannorein]

Glad it got resolved!