bench

rebound/simulation.py

@@ -457,6 +457,9 @@ def __repr__(self):
                 ("min_dt", c_double),
                 ("epsilon_global", c_uint),
                 ("dt_mode", c_uint),
+                ("rejected", c_long),
+                ("mindt", c_double),
+                ("safety", c_double),
                 ("_iterations_max_exceeded", c_ulong),
                 ("_allocatedN", c_int),
                 ("_at", POINTER(c_double)),

src/integrator_ias15.c

@@ -601,6 +601,7 @@ static int reb_integrator_ias15_step(struct reb_simulation* r) {
         if (fabs(dt_new)<r->ri_ias15.min_dt) dt_new = copysign(r->ri_ias15.min_dt,dt_new);
 
         if (fabs(dt_new/dt_done) < safety_factor) { // New timestep is significantly smaller.
+            r->ri_ias15.rejected += 1;
             // Reset particles
             for(int k=0;k<N;++k) {
                 int mk = map[k];
@@ -625,9 +626,14 @@ static int reb_integrator_ias15_step(struct reb_simulation* r) {
             return 0; // Step rejected. Do again. 
         }       
         if (fabs(dt_new/dt_done) > 1.0) {   // New timestep is larger.
-            if (dt_new/dt_done > 1./safety_factor) dt_new = dt_done /safety_factor; // Don't increase the timestep by too much compared to the last one.
+            if (dt_new/dt_done > 1./safety_factor){
+                r->ri_ias15.rejected += 1;
+                dt_new = dt_done /safety_factor; // Don't increase the timestep by too much compared to the last one.
+            }
         }
         r->dt = dt_new;
+        #define MIN(a, b) ((a) > (b) ? (b) : (a))
+        r->ri_ias15.mindt = MIN(r->ri_ias15.mindt, r->dt);
     }
 
     // Find new position and velocity values at end of the sequence

src/rebound.c

@@ -607,6 +607,9 @@ void reb_init_simulation(struct reb_simulation* r){
     r->ri_ias15.epsilon_global  = 1;
     r->ri_ias15.iterations_max_exceeded = 0;    
     r->ri_ias15.dt_mode = 1;
+    r->ri_ias15.rejected = 0;
+    r->ri_ias15.mindt = 1e300;
+    r->ri_ias15.safety = 0.25;
 
     // ********** SEI
     r->ri_sei.OMEGA     = 1;

src/rebound.h

@@ -120,6 +120,9 @@ struct reb_simulation_integrator_ias15 {
     double min_dt;
     unsigned int epsilon_global;
     unsigned int dt_mode;
+    long rejected;
+    double mindt;
+    double safety;
 
     // Internal use
     unsigned long iterations_max_exceeded; // Counter how many times the iteration did not converge.