MPI allocate only 1 particle in my simulation - Impossible to join multiple snapshots into one in MPI simulations

[Francyrad]

Dear @hannorein thank you in the first place for updating mpi. The examples work very well on ARM64

I wanted to implement mpi in my simulation and following the examples and guidance. However, when I run my simulation, it turns out I only have one particle. My code has some problems with particle distribution at various mpi nodes, but I can’t find the problem.

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include "rebound.h"
#include "tools.h"
#include "output.h"

void heartbeat(struct reb_simulation* const r);

const double G = 6.67430e-11;
const double M_earth = 5.972e24;
const double R_earth = 6371000.0;
const double rho_particle = 3300.0; // kg/m^3
const double r_inner = 6372000.0;
const double r_outer = 25484000.0;
const int N_particles = 10000;

const double M_moon = 7.342e22; // kg (massa lunare)
const double total_mass = M_moon; // massa totale delle particelle = 1 massa lunare

double coefficient_of_restitution_constant(const struct reb_simulation* const r, double v) {
    // v is the normal impact velocity.
    // Here, we just use a constant coefficient of restitution
    return 0.6;
}

int main(int argc, char* argv[]){
    struct reb_simulation* const r = reb_create_simulation();
    // Setup constants
    r->integrator    = REB_INTEGRATOR_LEAPFROG;
    r->gravity    = REB_GRAVITY_TREE;
    r->boundary    = REB_BOUNDARY_OPEN;
    r->collision = REB_COLLISION_TREE;
    r->coefficient_of_restitution = coefficient_of_restitution_constant;
    r->collision_resolve = reb_collision_resolve_hardsphere;
    r->opening_angle2    = 1.5;        // This constant determines the accuracy of the tree code gravity estimate.
    r->G         = G;
    r->softening     = 0.02;        // Gravitational softening length
    r->dt         = 3e-2;        // Timestep
    const double boxsize = 10.2;
    // Setup root boxes for gravity tree.
    // Here, we use 2x2=4 root boxes (each with length 'boxsize')
    // This allows you to use up to 8 MPI nodes.
    reb_configure_box(r,boxsize,4,2,1);

    // Initialize MPI
    // This can only be done after reb_configure_box.
    reb_mpi_init(r);

    // Setup particles only on master node
    double M_particle = total_mass / N_particles;
    double R_particle = pow((3.0 * M_particle) / (4.0 * M_PI * rho_particle), 1.0 / 3.0) / 1000.0; // Raggio in km

    struct reb_particle Earth = {0};
    Earth.m = M_earth;
    Earth.r = R_earth;
    if (r->mpi_id==0){
        reb_add(r, Earth);
        for (int i=0;i<N_particles;i++){
            struct reb_particle pt = {0};
            double a    = reb_random_powerlaw(r, boxsize/10.,boxsize/2./1.2,-1.5);
            double phi     = reb_random_uniform(r, 0,2.*M_PI);
            pt.x         = a*cos(phi);
            pt.y         = a*sin(phi);
            pt.z         = a*reb_random_normal(r, 0.001);
            double mu     = Earth.m + M_particle * (pow(a,-3./2.)-pow(boxsize/10.,-3./2.))/(pow(boxsize/2./1.2,-3./2.)-pow(boxsize/10.,-3./2.));
            double vkep     = sqrt(r->G*mu/a);
            pt.vx         =  vkep * sin(phi);
            pt.vy         = -vkep * cos(phi);
            pt.vz         = 0;
            pt.m         = M_particle;
            reb_add(r, pt);
        }
    }
    r->heartbeat = heartbeat;

#ifdef OPENGL
    // Hack to artificially increase particle array.
    // This cannot be done once OpenGL is activated.
    r->allocatedN *=8;
    r->particles = realloc(r->particles,sizeof(struct reb_particle)*r->allocatedN);
#endif // OPENGL

    // Start the integration
    reb_integrate(r, INFINITY);

    // Cleanup
    reb_mpi_finalize(r);
    reb_free_simulation(r);
    reb_simulationarchive_automate_interval(r, "archive.bin", 10.);
}

void heartbeat(struct reb_simulation* const r){
    if (reb_output_check(r,10.0*r->dt)){
        reb_output_timing(r,0);
    }
}

I edited the accretion disk example. Where i'm wrong?

Thank you in advance for the help Francesco

[hannorein]

Two issues I can see:

• Your earth radius is significantly larger than the simulation box. That doesn't make sense.
• You cannot call reb_simulationarchive_automate_interval after you free the simulation. That will lead to a segmentation fault.

[hannorein]

One more

• Also, your timestep is way too large. Your length scales are of the order of unit, your velocities are of the order of 1e7, but you have a timestep of 3e-2. It needs to be at least 1e-8 or 1e-9. In short, I think you need to revisit the units you're using (see https://rebound.readthedocs.io/en/latest/units/ if this is confusing).

[Francyrad]

Thank you professor, now it works. The only problem is that when i create my snapshots they are divided by mpi for obvious reasons. Is there not any way to avoid this? it seems impossible to put the archive function after mpi command.

[hannorein]

Glad it works. Yes, there will be one archive per node. Just create a separate file name for each node (using sprintf) and pass that to automate_simulationarchive then you'll get one file per node. Then use the same filenames when loading the data back in.

[Francyrad]

It already creates archivebin* files. My question now, when i have to create my pictures, i have to load the 4 archives and then arrange them merge them in some way?

[hannorein]

If you just want to visualize them and memory is not an issue, then I would simply load them into one simulation. You could do something like this:

struct reb_simualtion* main = reb_create_simulation();
for loop over all files{
   struct reb_simulation* tmp = reb_create_simulation_from_simulationarchive(archive, -1); 
   for (int i=0; i<tmp->N; i++){
     reb_add(main, tmp->particles[i]);
   }
}

[Francyrad]

I'm sorry for my late response. This is what i did:

void combine_snapshots(const int num_files) {
    char filename[256];
    char combined_filename[256];

    // Creiamo una nuova simulazione principale
    struct reb_simulation* main = reb_create_simulation();

    // Carichiamo ogni snapshot nell'archivio principale
    for (int file_num = 0; file_num < num_files; file_num++) {
        sprintf(filename, "archive%d.bin", file_num);
        struct reb_simulation* tmp = reb_create_simulation_from_simulationarchive(filename, -1);

        for (int i = 0; i < tmp->N; i++) {
            reb_add(main, tmp->particles[i]);
        }

        reb_free_simulation(tmp);
    }

    // Salva la simulazione combinata in un unico file binario
    sprintf(combined_filename, "combined_archive.bin");
    reb_output_binary(main, combined_filename);

    // Libera la memoria della simulazione principale
    reb_free_simulation(main);
}

And then at the end:

        // Start the integration
        char filename[256]; // This will hold the filename. Adjust the size as needed.
        sprintf(filename, "archive.bin", r->mpi_id); // This will create a filename like "archive0.bin", "archive1.bin", etc.
        reb_simulationarchive_automate_interval(r, "archive.bin", 100.);
        reb_integrate(r, INFINITY);

        // Cleanup and create snapshot.

        reb_mpi_finalize(r);

        reb_free_simulation(r);

        combine_snapshots(8);  // Assuming you have 8 snapshots as per your description

        return 0;

    }

The simulation compiles, but when i run it, the live visualization appears and i can see 1/4 of my accretion disk, then it crashes, so i don't know if the solution works...

I'm i missing something? thank you for your work. I ask also if it is possible to add an example about combining snapshots using mpi

Thank you for your availability and work

[hannorein]

Do I understand it correctly that you're trying to use MPI and OpenGL at the same time? I would run the simulations with MPI. Then visualize (but not integrate) them later with MPI turned off.

[Francyrad]

Here I am again!

I've been able to join all the particles in my simulation, however, i can't see my central particle. I used this python script:

import rebound
import glob

def main():
    # Gather all archive files
    archive_files = sorted(glob.glob("archive.bin_*"))

    if not archive_files:
        print("No archive files found!")
        return

    # Create a new simulation
    main_sim = rebound.Simulation()

    # Try finding the central particle in any of the archive files
    central_particle_found = False

    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        for snapshot in sa:
            if len(snapshot.particles) > 0:
                main_sim.add(snapshot.particles[0])  # Adding only the central particle
                central_particle_found = True
                break
        if central_particle_found:
            break

    if not central_particle_found:
        print("No central particle found in any snapshot of any archive file!")
        return

    # Now loop through all files and add particles, but skip the central particle
    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        tmp_sim = sa[-1]

        for particle in tmp_sim.particles[1:]:  # Start from 1 to exclude the central body which we've already added
            main_sim.add(particle)

    # At this point, main_sim has the central body and all the particles from the last snapshot of each archive
    print(f"Number of combined particles: {len(main_sim.particles)}")

    # Save the combined simulation to a new archive
    main_sim.save("archive_combined.bin")

if __name__ == "__main__":
    main()

and then:

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

def main(): archive_file = "archive_combined.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 particle in sim.particles:
    radius = particle.r
    if particle.m > 1e20:  # If the particle's mass is greater than 1e20 kg
        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()

Where does the problem could came from? I don't understand why my central particles is not white with a mass larger than 1e20

Francesco

[hannorein]

The central particle will not remain at index 0. So your lines:

main_sim.add(snapshot.particles[0]) 

and

for particle in tmp_sim.particles[1:]: 

will not have the intended effect. Just add all the particles from each archive once.

[Francyrad]

thank you, it worked, but my central particle is still gray...

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

def main():
    archive_file = "archive_combined.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()

the script worked when the code was run with OpenMP

[hannorein]

The index of your central particle is not 0, so this statement doesn't have the intended effect:

  if i == 0:  # Assuming the first particle is the planet
            color = 'white'

You could check the radius or mass instead...

[Francyrad]

I'm still sorry for the disturb, but when my archive is been created and divided in 8 "archive.bin_0" etc, these archives are constantly updated and there is no accumulation. So I have only 1 snapshot at a time. I didn't find any solution to that, in the old code without mpi it worked:

#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#include <time.h>
#include <sys/time.h>
#include <sys/stat.h>
#include <string.h>
#include <omp.h>
#include "rebound.h"
#include "tools.h"
#include "output.h"

void heartbeat(struct reb_simulation* const r);
void heartbeat(struct reb_simulation* const r){
    if (reb_output_check(r,10.0*r->dt)){
        reb_output_timing(r,0);
    }
}

const double G = 6.67430e-11;
const double M_earth = 5.972e24;
const double R_earth = 6371000.0;
const double rho_particle = 3300.0; // kg/m^3
const double r_inner = 6372000.0;
const double r_outer = 25484000.0;
const double boxsize =  3 * r_outer;
const int N_particles = 10000;

const double M_moon = 7.342e22; // kg (massa lunare)
const double total_mass = M_moon; // massa totale delle particelle = 1 massa lunare
double random_uniform(double min, double max) {
    return min + (max - min) * ((double) rand() / (double) RAND_MAX);
}
double coefficient_of_restitution_constant(const struct reb_simulation* const r, double v) {
    // v is the normal impact velocity.
    // Here, we just use a constant coefficient of restitution
    return 0.6;
}

int main(int argc, char* argv[]){
    char* filename = "archive.bin";
    struct reb_simulation* r = NULL;

    // Try to load a simulation from the archive
    struct reb_simulationarchive* sa = reb_open_simulationarchive(filename);
    if (sa != NULL) {
        r = reb_create_simulation_from_simulationarchive(sa, -1);
        reb_close_simulationarchive(sa);
    }

    if (r == NULL) {
        printf("No simulation archive found. Creating new simulation.\n");
        r = reb_create_simulation();
    // Setup constants
    r->integrator    = REB_INTEGRATOR_LEAPFROG;
    r->gravity    = REB_GRAVITY_TREE;
    r->boundary    = REB_BOUNDARY_OPEN;
    r->collision = REB_COLLISION_TREE;
    r->coefficient_of_restitution = coefficient_of_restitution_constant;
    r->collision_resolve = reb_collision_resolve_hardsphere;
    r->opening_angle2    = 1.5;        // This constant determines the accuracy of the tree code gravity estimate.
    r->G         = G;
    r->softening     = 0.02;        // Gravitational softening length
    r->dt         = 1;        // Timestep

    // Setup root boxes for gravity tree.
    // Here, we use 2x2=4 root boxes (each with length 'boxsize')
    // This allows you to use up to 8 MPI nodes.
    reb_configure_box(r, 3* boxsize, 4, 2, 1);

    // Initialize MPI
    // This can only be done after reb_configure_box.
    reb_mpi_init(r);

    // Setup particles only on master node

    struct reb_particle Earth = {0};
    Earth.m = M_earth;
    Earth.r = R_earth;

    double M_particle = total_mass / N_particles;
    double R_particle = pow((3.0 * M_particle) / (4.0 * M_PI * rho_particle), 1.0 / 3.0) / 1000.0; // Raggio in km

    if (r->mpi_id==0){
        reb_add(r, Earth);
        for (int i=0;i<N_particles;i++){
            double radius = random_uniform(r_inner, r_outer);
            double theta = random_uniform(0, 2 * M_PI);

            struct reb_particle p = {0};
            p.m = M_particle;
            p.r = R_particle * 1000.0; // Raggio in metri
            p.x = radius * cos(theta);
            p.y = radius * sin(theta);
            p.z = 0;
            double v_kep = sqrt(G * M_earth / radius);
            p.vx = -v_kep * sin(theta);
            p.vy = v_kep * cos(theta);
            p.vz = 0;

            reb_add(r, p);
                        }
                    }
                } else {
                    printf("Found simulation archive. Loaded snapshot at t=%.16f.\n", r->t);
                }

    r->heartbeat = heartbeat;
    reb_simulationarchive_automate_interval(r, "archive.bin", 100.);
    reb_integrate(r, INFINITY);

    reb_mpi_finalize(r);
    reb_free_simulation(r);

    }

Still sorry for the disturb and thank you for the availability

[hannorein]

You're only looking at the last snapshot?

# Now loop through all files and add particles, but skip the central particle
    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        tmp_sim = sa[-1]
...

[Francyrad]

No :(, it just save the last snapshot and substitute the old one. So i have only that instead of having a super massive archive divided in 8 tha accumulate memory:

The files sizes doesn't increase, it's been replaced time by time so i have only one snapshot to read. I'm sure that my python script for reading files works because I can read multiple snapshots in the non-mpi version of my code. I wouldn't like this is a mpi bug.

Can you please run the code and check? here is everything you need: bug.zip

I'm infinitely in debt

[hannorein]

See if 691e6dda fixes your problem.

[Francyrad]

I substituted archive.c but nothing... It just replace the old files as usual

[hannorein]

Next try. See if the latest version on the main branch helps.

[Francyrad]

Thank you, now they accumulate, so the snapshots are save, however, when I try to combine them, i got only one (the first), so 1 mb of combined file, using this script:

import rebound
import glob

def main():
    # Gather all archive files
    archive_files = sorted(glob.glob("archive.bin_*"))

    if not archive_files:
        print("No archive files found!")
        return

    # Create a new simulation
    main_sim = rebound.Simulation()

    # Try finding the central particle in any of the archive files
    central_particle_found = False

    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        for snapshot in sa:
            if len(snapshot.particles) > 0:
                main_sim.add(snapshot.particles[0])  # Adding only the central particle
                central_particle_found = True
                break
        if central_particle_found:
            break

    if not central_particle_found:
        print("No central particle found in any snapshot of any archive file!")
        return

    # Now loop through all files and add particles, but skip the central particle
    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        tmp_sim = sa[-1]

        for particle in tmp_sim.particles[0:]:  # Start from 0 to include the central particle
            main_sim.add(particle)

    # At this point, main_sim has the central body and all the particles from the last snapshot of each archive
    print(f"Number of combined particles: {len(main_sim.particles)}")

    # Save the combined simulation to a new archive
    main_sim.save("archive_combined.bin")

if __name__ == "__main__":
    main()

Then i got this message and I can create the combined bin only for the first snapshot:

python3.11 combine.py
/opt/homebrew/lib/python3.11/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)

Error! Cannot add two particles with the same coordinates to the tree.

Error! Cannot add two particles with the same coordinates to the tree. Number of combined particles: 10002

[hannorein]

Oh well, then this has to wait for a little. I've been working on some changes on how the input/output is handled on the show_settings branch for a while. This you can try it out if you want to. Once these changes are merged into the main branch, I can look at your problem again...

[Francyrad]

Here i am! I downloaded the latest version, the bug is still present. MPI create 8 archives for my case, when I try to join them with my script, i get:

francyrad@MacBook-Pro-di-Francesco moon_sim % python3.11 combine.py
/opt/homebrew/lib/python3.11/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)

Error! Cannot add two particles with the same coordinates to the tree.

Error! Cannot add two particles with the same coordinates to the tree.
Number of combined particles: 10002
francyrad@MacBook-Pro-di-Francesco moon_sim % python3.11 moon_plot.py
Number of snapshots: 1
Processing snapshot 0
Snapshot 0 processed

It just read and join the first snapshot:

The point is that I didn't understand if this is a rebound problem or a problem of my script. I attach here the .zip with all the files:

bug_mpi.zip

Thank you for the help and support

[hannorein]

If I use the latest REBOUND version and run this script, I seem to be getting the right result. Am I missing something?

import rebound
import glob

def main():
    # Gather all archive files
    archive_files = sorted(glob.glob("archive.bin_*"))

    if not archive_files:
        print("No archive files found!")
        return

    # Create a new simulation
    main_sim = rebound.Simulation()

    # Now loop through all files and add particles
    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        if sa[-1].N>0:
            for particle in sa[-1].particles:  
                main_sim.add(particle)

    # At this point, main_sim has the central body and all the particles from the last snapshot of each archive
    print(f"Number of combined particles: {main_sim.N}")

    # Save the combined simulation to a new archive
    main_sim.save("archive_combined.bin")

if __name__ == "__main__":
    main()

[Francyrad]

If I use the latest REBOUND version and run this script, I seem to be getting the right result. Am I missing something?

import rebound
import glob

def main():
    # Gather all archive files
    archive_files = sorted(glob.glob("archive.bin_*"))

    if not archive_files:
        print("No archive files found!")
        return

    # Create a new simulation
    main_sim = rebound.Simulation()

    # Now loop through all files and add particles
    for archive_file in archive_files:
        sa = rebound.SimulationArchive(archive_file)
        if sa[-1].N>0:
            for particle in sa[-1].particles:  
                main_sim.add(particle)

    # At this point, main_sim has the central body and all the particles from the last snapshot of each archive
    print(f"Number of combined particles: {main_sim.N}")

    # Save the combined simulation to a new archive
    main_sim.save("archive_combined.bin")

if __name__ == "__main__":
    main()

When I run the above script the archove bins are combined, but i get only and only the first snapshot, while all the others are missing...

You can see from the size. The problem is that it doesn't combine all the snapshot, but only the first one.

[hannorein]

sa[-1] selects the last snapshot. You can change it to sa[0] for the first, sa[1] the second, etc...

[Francyrad]

Thank you for the intuition, I solved with this:

import rebound
import glob

def main():
    # Gather all archive files
    archive_files = sorted(glob.glob("archive.bin_*"))

    if not archive_files:
        print("No archive files found!")
        return

    # Path for the combined archive
    combined_archive_path = "archive_combined.bin"

    # Initialize a new simulation for creating the combined archive
    combined_sim = rebound.Simulation()

    # Set the simulation archive filename and the snapshot interval
    combined_sim.automateSimulationArchive(combined_archive_path, interval=1, deletefile=True)

    # Find the number of snapshots in the first archive file
    num_snapshots = len(rebound.SimulationArchive(archive_files[0]))

    # Iterate over all snapshots
    for i in range(num_snapshots):
        # Clear the particles in the combined_sim
        while combined_sim.N > 0:
            combined_sim.remove(0)

        # Iterate over all archive files and add the i-th snapshot from each file to combined_sim
        for archive_file in archive_files:
            sa = rebound.SimulationArchive(archive_file)
            snapshot = sa[i]
            if snapshot.N > 0:
                for particle in snapshot.particles:
                    combined_sim.add(particle)

        # Manually set the time of the snapshot and integrate the simulation to save the snapshot
        combined_sim.t = snapshot.t
        combined_sim.integrate(snapshot.t)

    print(f"Combined archive created at {combined_archive_path}")

if __name__ == "__main__":
    main()

Please add this scripts in the examples, in this way it will be possible for everyone to use mpi without other troubles. Thank you for the help!