An improved adaptive timestep for IAS15

changelog.md

@@ -4,6 +4,14 @@ This changelog only includes the most important changes in recent updates. For a
 
 ## Version 3.x
 
+### Version 3.27.0
+* In python, Simulation and Particle objects are now picklable. Just like loading Simulations from a binary file, function pointers will need to be re-set manually after unpickling.
+* The difference between simulations can now be printed out in a human readable form. Python syntax: `sim.diff(sim2)`. C syntax: `reb_diff_simulations(sim2, sim1, 1)`.
+* Reading SimulationArchives with version < 2 is no longer supported.
+* The POSIX function fmemopen() is now required to compile REBOUND. This should not affect many users. However, if you are using macOS, the version needs to be >= 10.13 (this version of macOS, High Sierra, was released in 2017). 
+* Internal changes on how SimulationArchives are written. 
+* Internal variable names that represent the size of allocated buffers now consistently include the name `allocated_N`.
+
 ### Version 3.26.3
 * A few more changes to reduce the number of compiler warnings. This should not affect any calculation.
 

docs/binaryformat.md

@@ -1,9 +1,11 @@
 # Binary Format
 
-REBOUND comes with its one binary format.
-The binary format allows you to store a current simulation state to a file.
-The Simulation Archive is an extension of that format which allows you to store multiple snapshots in one file.
+REBOUND comes with its own binary format.
+The binary format allows you to store a current simulation state to a file or to memory.
+The binary format is also used when you make a copy of a simulation or when you compare two similations with each other.
+The Simulation Archive is an extension of the binary format which allows you to store multiple snapshots of a simulation in one file.
 This page explains the details of the binary format.
+It is mainly intended for people who wish to extend the built-in REBOUND functionality.
 You do not need to know those details if you're only working with binary files to save and load simulations.
 
 REBOUND uses two structures for the binary files:
@@ -43,7 +45,7 @@ You create a binary file if you save a simulation
     // ... setup simulation ...
     sim.save("snapshot.bin")
     ```
-Such a binary file with one snapshot is simply a set of `reb_binaryfield`s followed by one `reb_simulationarchive_blob` at the end:
+Such a binary file with one snapshot is simply a set of `reb_binaryfield`s followed by one `reb_simulationarchive_blob` at the end, for example:
 
 ```
 reb_binary_field:
@@ -71,8 +73,11 @@ reb_simulationarchive_blob:
 ```
 
 Each of the binary fields provides the context (type and size) for the data that immediately follows the field.
-The type is an integer defined in the `rebound.h` file as `enum REB_BINARY_FIELD_TYPE`.
-The last binary field is of type `END` to indicate that the snapshot ends here. 
+The type is an integer defined in the `reb_binary_field_descriptor_list` (see below).
+The last binary field of type `9999` (`end`) to indicate that the snapshot ends here. 
+
+!!! note
+    Before version 3.27 data was encoded using the enum `REB_BINARY_FIELD_TYPE` instead of `reb_binary_field_descriptor_list`.
 
 
 ## SimulationArchive file (multiple snapshots)
@@ -143,3 +148,25 @@ reb_simulationarchive_blob:
 The offsets are also used as a sort of checksum to detect if a binary file has been corrupted (for example because a user ran out of disk space). 
 If a binary file is corrupted, REBOUND attempts some magic and will recover the last snapshot which does not appear corrupted.
 You will see a warning message when that happens and should proceed with caution (make a backup!). 
+
+
+## Binary Field Descriptor
+
+REBOUND maintains a list of fields it needs to input/output in order to restore a simulation. 
+This list is of type `struct reb_binary_field_descriptor[]` and defined in `output.c` as `reb_binary_field_descriptor_list`.
+A single struct `reb_binary_field_descriptor` contains the information to input/output one REBOUND field, for example the current simulation time `t`:
+
+```c
+    struct reb_binary_field_descriptor fd_t = { 0, REB_DOUBLE, "t", offsetof(struct reb_simulation, t), 0, 0};
+```
+The first number is a unique identifier (in this case 0). The second entry is the type of data, in this case a single double precision floating point number. The third entry is a string used to identify the field. This is only used when generating human-readable output and is typically the same as the variable name in C. The next entry is the offset of where this variable is stored relative to the beginning of the simulation structure. 
+
+REBOUND also supports array like fields. For example consider the `particles` field:
+```c
+    struct reb_binary_field_descriptor fd_particles = { 85, REB_POINTER, "particles", offsetof(struct reb_simulation, particles), offsetof(struct reb_simulation, N), sizeof(struct reb_particle)};
+```
+
+The second to last entry lists the offset of the a variable in the `reb_simulation` structure that determines the number of array elements. In this case the number of partiles. The last entry is the size of a single element. In this case, the size of one `reb_particle`.
+
+If you add an additional field to the `reb_simulation` struct and you want to write it to a binary file and read it back in, then you need to add an entry to `reb_binary_field_descriptor_list`.
+

examples/binary_field_type_reverse/Makefile

@@ -0,0 +1,22 @@
+export OPENGL=0
+include ../../src/Makefile.defs
+
+all: librebound
+	@echo ""
+	@echo "Compiling problem file ..."
+	$(CC) -I../../src/ -Wl,-rpath,./ $(OPT) $(PREDEF) problem.c -L. -lrebound $(LIB) -o rebound
+	@echo ""
+	@echo "REBOUND compiled successfully."
+
+librebound: 
+	@echo "Compiling shared library librebound.so ..."
+	$(MAKE) -C ../../src/
+	@-rm -f librebound.so
+	@ln -s ../../src/librebound.so .
+
+clean:
+	@echo "Cleaning up shared library librebound.so ..."
+	@-rm -f librebound.so
+	$(MAKE) -C ../../src/ clean
+	@echo "Cleaning up local directory ..."
+	@-rm -vf rebound

examples/binary_field_type_reverse/problem.c

@@ -0,0 +1,43 @@
+/**
+ * A very simple test problem
+ * 
+ * We first create a REBOUND simulation, then we add 
+ * two particles and integrate the system for 100 time 
+ * units.
+ */
+#include "rebound.h"
+#include "output.h"
+#include <stdio.h>
+#include <stdlib.h>
+
+int main(int argc, char* argv[]) {
+    struct reb_simulation* r = reb_create_simulation();
+
+    reb_add_fmt(r, "m", 1.);                // Central object
+    reb_add_fmt(r, "m a e", 1e-3, 1., 0.1); // Jupiter mass planet
+    reb_add_fmt(r, "a e", 1.4, 0.1);        // Massless test particle 
+
+    reb_steps(r,1);
+
+    char* buf = NULL;
+    size_t size = 0;
+    reb_output_binary_to_stream(r, &buf, &size);
+
+
+
+    struct reb_simulationarchive* sa = malloc(sizeof(struct reb_simulationarchive));
+    reb_read_simulationarchive_from_buffer_with_messages(sa, buf, size, NULL, NULL);
+
+
+    struct reb_simulation* r2 = reb_create_simulation();
+    enum reb_input_binary_messages w =0;
+    reb_create_simulation_from_simulationarchive_with_messages(r2,sa,-1,&w);
+
+    reb_close_simulationarchive(sa);
+    free(buf);
+    printf("sizeof(size_t)= %ld\n", sizeof(unsigned int));
+
+    reb_free_simulation(r);
+    reb_free_simulation(r2);
+}
+

examples/closeencounter/problem.c

@@ -15,7 +15,7 @@
 #include "rebound.h"
 
 void heartbeat(struct reb_simulation* r);
-
+double E0;
 int main(int argc, char* argv[]){
     struct reb_simulation* r = reb_create_simulation();
     r->dt         = 0.01*2.*M_PI;        // initial timestep
@@ -40,11 +40,14 @@ int main(int argc, char* argv[]){
         reb_add(r, planet); 
     }
     reb_move_to_com(r);        // This makes sure the planetary systems stays within the computational domain and doesn't drift.
+    E0 = reb_tools_energy(r);
     reb_integrate(r, INFINITY);
 }
 
 void heartbeat(struct reb_simulation* r){
     if (reb_output_check(r, 10.*2.*M_PI)){  
         reb_output_timing(r, 0);
     }
+    double E1 = reb_tools_energy(r);
+    printf("%e\n", fabs((E1-E0)/E0));
 }

examples/mpi_unittests/problem.c

@@ -60,9 +60,50 @@ void test_twobody(){
     assert(fabs(o.a-1.)<1e-3);
     assert(fabs(o.e-0.1)<1e-2);
 
+    printf("Checking input/output...\n");
+
+    com = reb_get_com(r); // Need to call this on all machines. 
+    if (r->mpi_id==0){
+        for (int i=0; i<10; i++){
+            reb_add_fmt(r, "m a primary hash", 0.01, 2.0+0.1*i, com, i);
+        }
+    }
+    reb_steps(r, 1);
+    {
+        // Delete any previous files
+        char filename[1024];
+        sprintf(filename, "out.bin_%d", r->mpi_id);
+        remove(filename);
+    }
+    reb_simulationarchive_snapshot(r, "out.bin");
+    reb_steps(r, 1);
+    reb_simulationarchive_snapshot(r, "out.bin");
+    reb_steps(r, 10);
+
+    struct reb_simulationarchive* sa = reb_open_simulationarchive("out.bin");
+    assert(sa->nblobs==2);
+    struct reb_simulation* r2 = reb_create_simulation_from_simulationarchive(sa,-1);
+    reb_steps(r2, 10);
+
+    assert(r->N == r2->N);
+    assert(r->t == r2->t);
+
+    // Order of particles will be different. Need to compare them by hash
+    for(int i=0; i<10; i++){
+        struct reb_particle p1 = reb_get_remote_particle_by_hash(r, i);
+        struct reb_particle p2 = reb_get_remote_particle_by_hash(r2, i);
+        assert(p1.x==p2.x);
+        assert(p1.y==p2.y);
+        assert(p1.z==p2.z);
+    }
+
+
+
+
     printf("Cleanup...\n");
     reb_mpi_finalize(r);
     reb_free_simulation(r); 
+    reb_free_simulation(r2); 
 
 
 }

examples/selfgravity_disc_mpi/problem.c

@@ -76,8 +76,8 @@ int main(int argc, char* argv[]){
 #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);
+    r->allocated_N *=8;
+    r->particles = realloc(r->particles,sizeof(struct reb_particle)*r->allocated_N);
 #endif // OPENGL
 
     // Start the integration

examples/shearing_sheet_mpi/problem.c

@@ -88,8 +88,8 @@ int main(int argc, char* argv[]) {
 #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);
+    r->allocated_N *=8;
+    r->particles = realloc(r->particles,sizeof(struct reb_particle)*r->allocated_N);
 #endif // OPENGL
 
     // Start the integration

rebound/__init__.py

@@ -84,7 +84,7 @@ class ParticleNotFound(Exception):
     pass
 
 from .tools import hash, mod2pi, M_to_f, E_to_f, M_to_E, spherical_to_xyz, xyz_to_spherical
-from .simulation import Simulation, Orbit, Variation, reb_simulation_integrator_saba, reb_simulation_integrator_whfast, reb_simulation_integrator_sei, reb_simulation_integrator_mercurius, reb_simulation_integrator_ias15, ODE, Rotation, Vec3d, _Vec3d
+from .simulation import Simulation, Orbit, Variation, reb_simulation_integrator_saba, reb_simulation_integrator_whfast, reb_simulation_integrator_sei, reb_simulation_integrator_mercurius, reb_simulation_integrator_ias15, ODE, Rotation, Vec3d, _Vec3d, binary_field_descriptor_list
 from .particle import Particle
 from .plotting import OrbitPlot, OrbitPlotSet
 from .simulationarchive import SimulationArchive
@@ -97,4 +97,4 @@ def __init__(self, processes=None, initializer=None, initargs=(), **kwargs):
 else:
     from .interruptible_pool import InterruptiblePool
 
-__all__ = ["__libpath__", "__version__", "__build__", "__githash__", "SimulationArchive", "Simulation", "Orbit", "OrbitPlot", "OrbitPlotSet", "Particle", "SimulationError", "Encounter", "Collision", "Escape", "NoParticles", "ParticleNotFound", "InterruptiblePool","Variation", "reb_simulation_integrator_whfast", "reb_simulation_integrator_ias15", "reb_simulation_integrator_saba", "reb_simulation_integrator_sei","reb_simulation_integrator_mercurius", "clibrebound", "mod2pi", "M_to_f", "E_to_f", "M_to_E", "ODE", "Rotation", "Vec3d", "spherical_to_xyz", "xyz_to_spherical"]
+__all__ = ["__libpath__", "__version__", "__build__", "__githash__", "SimulationArchive", "Simulation", "Orbit", "OrbitPlot", "OrbitPlotSet", "Particle", "SimulationError", "Encounter", "Collision", "Escape", "NoParticles", "ParticleNotFound", "InterruptiblePool","Variation", "reb_simulation_integrator_whfast", "reb_simulation_integrator_ias15", "reb_simulation_integrator_saba", "reb_simulation_integrator_sei","reb_simulation_integrator_mercurius", "clibrebound", "mod2pi", "M_to_f", "E_to_f", "M_to_E", "ODE", "Rotation", "Vec3d", "spherical_to_xyz", "xyz_to_spherical","binary_field_descriptor_list"]

rebound/particle.py

@@ -1,4 +1,4 @@
-from ctypes import Structure, c_double, c_int, byref, memmove, sizeof, c_uint32, c_uint, c_ulong
+from ctypes import Structure, c_double, c_int, byref, memmove, sizeof, c_uint32, c_uint, c_ulong, c_char_p, string_at, POINTER, c_char
 from . import clibrebound, E_to_f, M_to_f, mod2pi
 import math
 import ctypes.util
@@ -156,6 +156,25 @@ def __init__(self, simulation=None, particle=None, m=None, x=None, y=None, z=Non
 
         """        
 
+        # Unpickling
+        binarydata = None
+        if isinstance(simulation, bytes):
+            # simulation is actually a bytes array with the particle data
+            binarydata = simulation
+        if isinstance(particle, bytes):
+            # simulation is actually a bytes array with the particle data
+            binarydata = particle
+        if binarydata is not None:
+            if len(binarydata) != sizeof(self):
+                raise ValueError("Binary particle data does not have the right size.")
+            buft = c_char * len(binarydata)
+            buf = buft.from_buffer_copy(binarydata)
+            memmove(byref(self), byref(buf), sizeof(self))
+            self.c = 0
+            self.sim = 0
+            self.ap = 0
+            return
+
         if Omega == "uniform":
             Omega = random.vonmisesvariate(0.,0.) 
         if omega == "uniform":
@@ -401,6 +420,10 @@ def copy(self):
         memmove(byref(np), byref(self), sizeof(self))
         return np
 
+# Pickling method
+    def __reduce__(self):
+        return (Particle, (string_at(byref(self), size=sizeof(self)),))
+
     def calculate_orbit(self, primary=None, G=None):
         """ 
         Returns a rebound.Orbit object with the keplerian orbital elements
@@ -572,6 +595,14 @@ def sample_orbit(self, Npts=100, primary=None, samplingAngle=None, duplicateEndp
 
     # Simple operators for particles.
 
+    def __eq__(self, other):
+        # This ignores the pointer values
+        if not isinstance(other,Particle):
+            return NotImplemented
+        clibrebound.reb_diff_particles.restype = c_int
+        ret = clibrebound.reb_diff_particles(self, other)
+        return not ret
+
     def __pow__(self, other):
         if not isinstance(other, Particle):
             return NotImplemented