Updates to SeBa (#1128)

* added two tests to test_seba.py to check earlier seba interface problems. And relaired test_seba: should now work as expected.

* assure that new bug-corrected version of SeBa is downloaded (see SeBa issue a6f4), and make SeBa interface print zero-age binary parameters.

* Update download.py

* Update download.py

---------

Co-authored-by: Simon Portegies Zwart <spz@strw.leidenuniv.nl>
Co-authored-by: Steven Rieder <steven+github@rieder.nl>

Author: 3 people

src/amuse/community/seba/download.py

@@ -93,7 +93,7 @@ def new_option_parser():
     result = OptionParser()
     result.add_option(
         "--seba-version",
-        default='2d8088ad03a4323514780e19e5895fbcac42e0ec',
+        default='a6f4b64f852249a0b66d52aa6f276be517dccb64',
         dest="seba_version",
         help="SeBa commit hash to download",
         type="string"

src/amuse/community/seba/interface.cc

@@ -42,8 +42,8 @@ local void addbinary(
 
             id = bi->get_index();
 
-            // cerr << "Adding binary to "<< id << " at time = "
-            //      << stellar_time << endl;
+            cerr << "Adding binary to "<< id << " at time = "
+                 << stellar_time << endl;
 
             double_star* new_double
             = new_double_star(bi, sma, ecc, stellar_time, id, type);
@@ -53,7 +53,8 @@ local void addbinary(
             // Give the new binary the old star_story.
 
             new_double->set_star_story(old_story);
-
+	    // Assure that SeBa prints zero-age binary parameters (SPZ+FK Febr2025).
+	    new_double->dump("SeBa.data", true);
 
         }
         else {
@@ -458,6 +459,7 @@ int new_particle(int * index_of_the_star, double mass){
     }
 
     addstar(new_node, seba_time, start_type, seba_metallicity, 0, false, start_type, mass, mass, 0, 0, 0);
+
     new_node->get_starbase()->set_time_offset(seba_time);
     *index_of_the_star = next_seba_id;
 
@@ -919,6 +921,15 @@ int new_binary(
     child2->set_elder_sister(child1);
     child1->set_parent(new_node);
     child2->set_parent(new_node);
+    // We cannot access the stellar parameters from here, because the starbase class is in between.
+    // SPZ*FK Febr.2025
+    //cerr << "Set stelar id in new binary interface."<<endl;
+    //cerr << "Current interface "<< child1->get_star_id()<<"and "<< child1->get_star_id()<<endl;
+    //child2->set_star_id(1);
+    // set stellar id's
+    //child1->set_star_id(0);
+    //child2->set_star_id(1);
+    
 
     addbinary(new_node, seba_time, binary_start_type, semi_major_axis, eccentricity);
 

src/amuse/test/suite/codes_tests/test_seba.py

@@ -134,8 +134,8 @@ def test5(self):
         print(mass)
         self.assertEqual(error, 0)
         self.assertAlmostRelativeEqual(mass[0], 1.0, 6)
-        self.assertAlmostRelativeEqual(mass[1], 0.62973, 4)
-        self.assertAlmostRelativeEqual(mass[2], 0.75012, 4)
+        self.assertAlmostRelativeEqual(mass[1], 0.6306, 4)
+        self.assertAlmostRelativeEqual(mass[2], 0.7408, 4)
 
         instance.stop()
 
@@ -155,7 +155,7 @@ def test6(self):
         mass, error = instance.get_mass(index)
         print(mass)
         self.assertEqual(error, 0)
-        self.assertAlmostRelativeEqual(mass, [1.0, 0.62973, 0.75072], 4)
+        self.assertAlmostRelativeEqual(mass, [1.0, 0.6306, 0.7408], 4)
 
         instance.stop()
 
@@ -284,7 +284,7 @@ def test9(self):
         self.assertEqual(error, 0)
         mass, error = instance.get_mass(1)
         self.assertEqual(error, 0)
-        self.assertAlmostRelativeEqual(mass, 2.98777, 4)
+        self.assertAlmostRelativeEqual(mass,  2.9887, 4)
         mass, error = instance.get_mass(2)
         self.assertEqual(error, 0)
         self.assertAlmostRelativeEqual(mass, 0.29999, 4)
@@ -293,7 +293,7 @@ def test9(self):
         self.assertEqual(error, 0)
         mass, error = instance.get_mass(1)
         self.assertEqual(error, 0)
-        self.assertAlmostRelativeEqual(mass, 0.902743, 4)
+        self.assertAlmostRelativeEqual(mass, 0.9019, 4)
         mass, error = instance.get_mass(2)
         self.assertEqual(error, 0)
         self.assertAlmostRelativeEqual(mass, 0.3, 4)
@@ -420,7 +420,7 @@ def test_evolution_of_close_binary_system(self):
             3.0000 | units.MSun,
             3.0000 | units.MSun,
             2.9983 | units.MSun,
-            2.9741 | units.MSun,
+            2.9797 | units.MSun,
             0.6710 | units.MSun,
             0.6596 | units.MSun,
         )
@@ -683,3 +683,53 @@ def test_supernova_stopping_condition_with_multiple_stars_of_equal_mass(self):
         self.assertAlmostRelativeEqual(instance.particles[0].mass, 1.2507 | units.MSun, 4)
         self.assertAlmostRelativeEqual(instance.particles[1].mass, 1.2507 | units.MSun, 4)
         self.assertAlmostRelativeEqual(instance.particles[2].mass, 0.5 | units.MSun, 4)
+
+    def test_restart_for_different_stellar_type(self):
+        """ Test restart SeBa for different stellar types. """
+        instance = self.new_instance_of_an_optional_code(SeBa)
+
+        #a very specific case, which went wrong in Torch
+        stars = Particles(2)
+        stars.mass = [69.08994562, 69.40049406] | units.MSun
+        stars.relative_mass = [69.44820097, 69.40049406] | units.MSun
+        stars.age = [3.2865, 3.2865] | units.Myr
+        stars.relative_age = [3.68045615, 3.63104588] | units.Myr
+        stars.stellar_type = [2, 1] | units.stellar_type
+        stars.core_mass = [29.90683271, 0.] | units.MSun
+        stars.radius = [944.70734099, 94.874201] | units.RSun
+        stars.luminosity = [1517964.62534381, 1468414.68887237] | units.LSun
+        
+        instance.particles.add_particles(stars)
+        instance.evolve_model(0.1 | units.day)
+        
+        self.assertAlmostRelativeEqual(instance.particles[0].age, 0 | units.Myr, 4)
+        self.assertAlmostRelativeEqual(instance.particles[1].age, 0 | units.Myr, 4)
+        self.assertAlmostRelativeEqual(instance.particles[0].relative_age,  3.68045615 | units.Myr, 4)
+        self.assertAlmostRelativeEqual(instance.particles[1].relative_age,  3.63104588 | units.Myr, 4)
+        self.assertAlmostRelativeEqual(instance.particles[0].radius,  944.70734099 | units.RSun, 4)
+        self.assertAlmostRelativeEqual(instance.particles[1].radius,  94.874201 | units.RSun, 4)
+        self.assertAlmostRelativeEqual(instance.particles[0].core_mass, 29.90683275 | units.MSun, 4)
+        self.assertAlmostRelativeEqual(instance.particles[1].core_mass, 0 | units.MSun, 4)
+        self.assertAlmostRelativeEqual(instance.particles[0].luminosity, 1517964.62534381 | units.LSun, 4)
+        self.assertAlmostRelativeEqual(instance.particles[1].luminosity,  1468414.68887237 | units.LSun, 4)
+
+
+    def test_adding_zero_age_star_and_recover_with_channel(self):
+        """ Test restart SeBa for different stellar types. """
+        instance = self.new_instance_of_an_optional_code(SeBa)
+
+        #a very specific case, which went wrong in Torch
+        
+        stars = Particles(1)
+        stars.mass = 20|units.MSun
+        stars.radius = 10|units.RSun
+
+        instance.particles.add_particle(stars[0])
+        instance.evolve_model(1|units.yr)
+
+        channel = stars.new_channel_to(instance.particles)
+        channel.copy()
+
+        self.assertAlmostRelativeEqual(instance.particles[0].mass, 19.9999999903 | units.MSun, 4)
+        self.assertAlmostRelativeEqual(instance.particles[0].radius,  5.99895518242 | units.RSun, 4)
+