Repository: hannorein/rebound (language: python)

Issue:
Additional forces double-counted in integrator_bs.c
Dear Hanno,

I believe I found an issue in the Burlisch-Stohr integrator in rebound. In particular, I think it is double-counting the accelerations from additional forces. 

I am sharing some code that illustrates this (https://gist.github.com/alekseygenerozov/fb28ca567a495100d0fce9b24170e8f4). This is a modified example problem from reboundx (in examples/central_force/problem.c). Basically the setup is such that (i) the additional force dominates (to better illustrate the problem) and (ii) the integration time is very short so you see the initial acceleration of the particle (the planet in this example). I am doing this test with the latest versions of rebound and reboundx

The acceleration in this example is almost exactly twice greater than expected. Here is the output I get 

t ax ay az: 1.000000e-15 3.125000e-04 4.784160e-25 0.000000e+00
x y z: 0.800000 0.000000 0.000000

The problem goes away when I change 

    sim->integrator=REB_INTEGRATOR_BS;

To 
    sim->integrator=REB_INTEGRATOR_IAS15;

I believe the issue is line 327 in src/integrator_bs.c. This line seems to be superfluous since additional_forces is also called in reb_update_acceleration. When I comment out this line the issue goes away. 

Best,

Aleksey



Discussion / hints:

Fix patch:
diff --git a/README.md b/README.md
index 84131057f..c5c9ea168 100644
--- a/README.md
+++ b/README.md
@@ -1,4 +1,4 @@
-[![Version](https://img.shields.io/badge/rebound-v3.19.8-green.svg?style=flat)](https://rebound.readthedocs.org)
+[![Version](https://img.shields.io/badge/rebound-v3.19.9-green.svg?style=flat)](https://rebound.readthedocs.org)
 [![PyPI](https://badge.fury.io/py/rebound.svg)](https://badge.fury.io/py/rebound)
 [![GPL](https://img.shields.io/badge/license-GPL-green.svg?style=flat)](https://github.com/hannorein/rebound/blob/main/LICENSE)
 [![Paper](https://img.shields.io/badge/arXiv-1110.4876-green.svg?style=flat)](https://arxiv.org/abs/1110.4876)
diff --git a/changelog.md b/changelog.md
index 3a89a0211..b8a0d145b 100644
--- a/changelog.md
+++ b/changelog.md
@@ -4,6 +4,9 @@ This changelog only includes the most important changes in recent updates. For a
 
 ## Version 3.x
 
+### Version 3.19.9
+* Two bugs fixed in the BS integrator. One was related to unitialized memory and the other to issues when the particle number changed.
+
 ### Version 3.19.5
 * Workaround for urllib support in pyodide added
 * Silent warning when InterruptiblePool is not available
diff --git a/setup.py b/setup.py
index 9a912d545..f77c83022 100644
--- a/setup.py
+++ b/setup.py
@@ -17,7 +17,7 @@
     ghash = subprocess.check_output(["git", "rev-parse", "HEAD"]).decode("ascii")
     ghash_arg = "-DGITHASH="+ghash.strip()
 except:
-    ghash_arg = "-DGITHASH=a6aeee6df105e0ee5216f55534f565727e5f7778" #GITHASHAUTOUPDATE
+    ghash_arg = "-DGITHASH=8219fda3fd538faa9ab7b4e5ffe23d46bbc3d244" #GITHASHAUTOUPDATE
 
 extra_link_args=[]
 if sys.platform == 'darwin':
@@ -64,7 +64,7 @@
     long_description = f.read()
 
 setup(name='rebound',
-    version='3.19.8',
+    version='3.19.9',
     description='An open-source multi-purpose N-body code',
     long_description=long_description,
     long_description_content_type="text/markdown",
diff --git a/src/integrator_bs.c b/src/integrator_bs.c
index fcce53a51..b50cb70e4 100644
--- a/src/integrator_bs.c
+++ b/src/integrator_bs.c
@@ -320,11 +320,6 @@ static void nbody_derivatives(struct reb_ode* ode, double* const yDot, const dou
         // Not needed for first step. Accelerations already calculated. Just need to copy them
         reb_integrator_bs_update_particles(r, y);
         reb_update_acceleration(r);
-        if (r->N_var){
-            reb_calculate_acceleration_var(r);
-        }
-        // Calculate non-gravity accelerations. 
-        if (r->additional_forces) r->additional_forces(r);
     }
 
     for (int i=0; i<r->N; i++){
diff --git a/src/rebound.c b/src/rebound.c
index 4033ed769..01938b9ba 100644
--- a/src/rebound.c
+++ b/src/rebound.c
@@ -65,7 +65,7 @@
 const int reb_max_messages_length = 1024;   // needs to be constant expression for array size
 const int reb_max_messages_N = 10;
 const char* reb_build_str = __DATE__ " " __TIME__;  // Date and time build string. 
-const char* reb_version_str = "3.19.8";         // **VERSIONLINE** This line gets updated automatically. Do not edit manually.
+const char* reb_version_str = "3.19.9";         // **VERSIONLINE** This line gets updated automatically. Do not edit manually.
 const char* reb_githash_str = STRINGIFY(GITHASH);             // This line gets updated automatically. Do not edit manually.
 
 static int reb_error_message_waiting(struct reb_simulation* const r);
diff --git a/version.txt b/version.txt
index afa2a6259..1db073dbe 100644
--- a/version.txt
+++ b/version.txt
@@ -1,1 +1,1 @@
-3.19.8
+3.19.9


Test patch:
diff --git a/rebound/tests/test_bs.py b/rebound/tests/test_bs.py
index f0e54367a..52b54feab 100644
--- a/rebound/tests/test_bs.py
+++ b/rebound/tests/test_bs.py
@@ -4,8 +4,32 @@
 import rebound.data
 import warnings
     
+def af(simp):
+    sim = simp.contents
+    x = sim.particles[0].x
+    y = sim.particles[0].y
+    z = sim.particles[0].z
+    r = math.sqrt(x*x+y*y+z*z)
+    sim.particles[0].ax -= x/(r*r*r)
+    sim.particles[0].ay -= y/(r*r*r)
+    sim.particles[0].az -= z/(r*r*r)
     
 class TestIntegratorBS(unittest.TestCase):
+    def test_bs_additionali_force_only(self):
+        sim = rebound.Simulation()
+        sim.additional_forces = af
+        sim.integrator = "bs"
+        eps = 1e-11
+        sim.ri_bs.eps_rel = eps
+        sim.ri_bs.eps_abs = eps
+        sim.add(m=0,x=1,vy=1)
+        sim.integrate(2.*math.pi)
+        self.assertLess(math.fabs(sim.particles[0].x-1.),5*eps)
+        self.assertLess(math.fabs(sim.particles[0].vy-1.),5*eps)
+        self.assertLess(math.fabs(sim.particles[0].y),5*eps)
+        self.assertLess(math.fabs(sim.particles[0].vx),5*eps)
+
+
     def test_bs_outersolarsystem(self):
         for eps in [1e-5, 1e-7, 1e-9, 1e-11]:
             sim = rebound.Simulation()