diff --git a/src/amuse/community/rebound/interface.cc b/src/amuse/community/rebound/interface.cc index 44d88cace9..ee9426bd20 100644 --- a/src/amuse/community/rebound/interface.cc +++ b/src/amuse/community/rebound/interface.cc @@ -161,8 +161,13 @@ int get_total_radius(double * radius){ return 0; } -int new_particle(int * index_of_the_particle, double mass, double x, - double y, double z, double vx, double vy, double vz, double radius, int code_index){ +int new_particle( + int * index_of_the_particle, + double mass, + double x, double y, double z, + double vx, double vy, double vz, + double radius, int code_index + ){ if(code_index < 0 || code_index >= (signed) codes.size()){ *index_of_the_particle=0; return -10; @@ -859,20 +864,20 @@ int _set_integrator(int value, int code_index){ code->integrator = reb_simulation::REB_INTEGRATOR_JANUS; break; case 9: - code->integrator = reb_simulation::REB_INTEGRATOR_WHFAST512; + code->integrator = reb_simulation::REB_INTEGRATOR_MERCURIUS; break; case 10: code->integrator = reb_simulation::REB_INTEGRATOR_SABA; break; case 11: - code->integrator = reb_simulation::REB_INTEGRATOR_MERCURIUS; - break; - case 12: code->integrator = reb_simulation::REB_INTEGRATOR_EOS; break; - case 13: + case 12: code->integrator = reb_simulation::REB_INTEGRATOR_BS; break; + case 21: + code->integrator = reb_simulation::REB_INTEGRATOR_WHFAST512; + break; default: code->integrator = reb_simulation::REB_INTEGRATOR_NONE; return -1; diff --git a/src/amuse/community/rebound/interface.py b/src/amuse/community/rebound/interface.py index 57c2b5016d..d167c6a796 100644 --- a/src/amuse/community/rebound/interface.py +++ b/src/amuse/community/rebound/interface.py @@ -195,11 +195,11 @@ def _get_integrator(): "whfast-helio": 6, # removed "none": 7, "janus": 8, - "whfast512": 9, + "mercurius": 9, "saba": 10, - "mercurius": 11, - "eos": 12, - "bs": 13, + "eos": 11, + "bs": 12, + "whfast512": 21, } def set_integrator(self, name, code_index=0): diff --git a/src/amuse/test/suite/codes_tests/test_rebound.py b/src/amuse/test/suite/codes_tests/test_rebound.py index da5933259d..695731cdca 100644 --- a/src/amuse/test/suite/codes_tests/test_rebound.py +++ b/src/amuse/test/suite/codes_tests/test_rebound.py @@ -1,24 +1,23 @@ -from amuse.community import * +import numpy as np +from amuse.datamodel import Particles +from amuse.units import nbody_system, units from amuse.test.amusetest import TestWithMPI - from amuse.community.rebound.interface import ReboundInterface from amuse.community.rebound.interface import Rebound -import math -try: - from matplotlib import pyplot - HAS_MATPLOTLIB = True -except ImportError: - HAS_MATPLOTLIB = False class ReboundInterfaceTests(TestWithMPI): - def test1(self): + def test_add_two_particles_serial_and_retrieve(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) instance.initialize_code() - res1, error = instance.new_particle(mass=11.0, radius=2.0, x=1.0, y=2.0, z=3.0, vx=4.0, vy=5.0, vz=6.0) - res2, error = instance.new_particle(mass=21.0, radius=5.0, x=10.0, y=0.0, z=0.0, vx=10.0, vy=0.0, vz=0.0) + res1, error = instance.new_particle( + mass=11.0, radius=2.0, x=1.0, y=2.0, z=3.0, vx=4.0, vy=5.0, vz=6.0 + ) + res2, error = instance.new_particle( + mass=21.0, radius=5.0, x=10.0, y=0.0, z=0.0, vx=10.0, vy=0.0, vz=0.0 + ) self.assertEqual(error, 0) self.assertEqual(res1, 0) @@ -46,11 +45,16 @@ def test1(self): instance.cleanup_code() instance.stop() - def test2(self): + def test_add_two_particles_parallel_and_retrieve(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) instance.initialize_code() - instance.new_particle([10, 20], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [0, 0], [1, 1]) + instance.new_particle( + [10, 20], + [0, 0], [0, 0], [0, 0], + [0, 0], [0, 0], [0, 0], + [1, 1] + ) retrieved_state = instance.get_state(0) self.assertEqual(10.0, retrieved_state['mass']) @@ -61,12 +65,16 @@ def test2(self): instance.cleanup_code() instance.stop() - def test3(self): - + def test_add_three_particles_and_retrieve(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) instance.initialize_code() - indices, error = instance.new_particle([10, 20, 30], [1, 2, 3], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [0, 0, 0], [1, 0, 1]) + indices, error = instance.new_particle( + [10, 20, 30], + [1, 2, 3], [0, 0, 0], [0, 0, 0], + [0, 0, 0], [0, 0, 0], [0, 0, 0], + [1, 0, 1] + ) n, error = instance.get_number_of_particles() self.assertEqual(n, 3) error = instance.delete_particle(indices[1]) @@ -82,39 +90,64 @@ def test3(self): instance.cleanup_code() instance.stop() - def test4(self): - + def test_set_integrators(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) instance.initialize_code() error = 0 - integrator = instance.get_integrator() - self.assertEqual(error, 0) - self.assertEqual("ias15", integrator) - error = instance.set_integrator("whfast") - self.assertEqual(error, 0) - integrator = instance.get_integrator() - self.assertEqual(error, 0) - self.assertEqual("whfast", integrator) + list_of_integrators = [ + "ias15", + "whfast", + "sei", + "leapfrog", + # "hermes", # removed + # "whfast-helio", # removed + "none", + "janus", + "mercurius", + "saba", + "eos", + "bs", + "whfast512", + ] + for i, integrator in enumerate(list_of_integrators): + print(i, integrator) + retrieved_integrator = instance.get_integrator() + self.assertEqual(error, 0) + self.assertEqual(integrator, retrieved_integrator) + try: + error = instance.set_integrator(list_of_integrators[i+1]) + self.assertEqual(error, 0) + except IndexError: + retrieved_integrator = instance.get_integrator() + self.assertEqual(error, 0) + self.assertEqual(integrator, retrieved_integrator) + instance.cleanup_code() instance.stop() - def test5(self): + def test_evolve_default_integrator(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) self.assertEqual(0, instance.initialize_code()) self.assertEqual(0, instance.commit_parameters()) # Set up an equal-mass binary on a circular orbit: - self.assertEqual([0, 0], list(instance.new_particle(0.5, 0.5, 0, 0, 0, 0.5, 0, 0.01).values())) - self.assertEqual([1, 0], list(instance.new_particle(0.5, -0.5, 0, 0, 0, -0.5, 0, 0.01).values())) + self.assertEqual( + [0, 0], + list(instance.new_particle(0.5, 0.5, 0, 0, 0, 0.5, 0, 0.01).values()) + ) + self.assertEqual( + [1, 0], + list(instance.new_particle(0.5, -0.5, 0, 0, 0, -0.5, 0, 0.01).values()) + ) self.assertEqual(0, instance.commit_particles()) - self.assertEqual(0, instance.evolve_model(math.pi)) + self.assertEqual(0, instance.evolve_model(np.pi)) for result, expected in zip(list(instance.get_position(0).values()), [-0.5, 0.0, 0.0, 0]): self.assertAlmostEqual(result, expected, 3) for result, expected in zip(list(instance.get_position(1).values()), [0.5, 0.0, 0.0, 0]): self.assertAlmostEqual(result, expected, 3) - self.assertEqual(0, instance.evolve_model(2 * math.pi)) + self.assertEqual(0, instance.evolve_model(2 * np.pi)) for result, expected in zip(list(instance.get_position(0).values()), [0.5, 0.0, 0.0, 0]): self.assertAlmostEqual(result, expected, 3) for result, expected in zip(list(instance.get_position(1).values()), [-0.5, 0.0, 0.0, 0]): @@ -123,18 +156,23 @@ def test5(self): self.assertEqual(0, instance.cleanup_code()) instance.stop() - def test6(self): + def test_fail_to_add_particle(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) self.assertEqual(0, instance.initialize_code()) self.assertEqual(0, instance.commit_parameters()) - index, err = instance.new_particle(0.5, 0.5, 0, 0, 0, 0.5, 0, 0.01, 1) + index, err = instance.new_particle( + 0.5, + 0.5, 0, 0, + 0, 0.5, 0, + 0.01, 1 + ) self.assertEqual(-10, err) instance.cleanup_code() instance.stop() - def test7(self): + def test_set_epsilon_squared(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) instance.initialize_code() @@ -147,7 +185,7 @@ def test7(self): instance.cleanup_code() instance.stop() - def test8(self): + def test_set_boundaries(self): instance = self.new_instance_of_an_optional_code(ReboundInterface) instance.initialize_code() @@ -170,23 +208,26 @@ def test8(self): class TestRebound(TestWithMPI): def new_system_of_sun_and_earth(self): - stars = datamodel.Stars(2) + stars = Particles(2) sun = stars[0] sun.mass = units.MSun(1.0) - sun.position = units.m(numpy.array((0.0, 0.0, 0.0))) - sun.velocity = units.ms(numpy.array((0.0, 0.0, 0.0))) + sun.position = units.m(np.array((0.0, 0.0, 0.0))) + sun.velocity = units.ms(np.array((0.0, 0.0, 0.0))) sun.radius = units.RSun(1.0) earth = stars[1] earth.mass = units.kg(5.9736e24) earth.radius = units.km(6371) - earth.position = units.km(numpy.array((149.5e6, 0.0, 0.0))) - earth.velocity = units.ms(numpy.array((0.0, 29800, 0.0))) + earth.position = units.km(np.array((149.5e6, 0.0, 0.0))) + earth.velocity = units.ms(np.array((0.0, 29800, 0.0))) return stars - def test1(self): - convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km) + def test_evolve_default_integrator(self): + convert_nbody = nbody_system.nbody_to_si( + 1.0 | units.MSun, + 149.5e6 | units.km, + ) interface = self.new_instance_of_an_optional_code(Rebound, convert_nbody) interface.initialize_code() @@ -224,6 +265,7 @@ def test1(self): interface.stop() def test2(self): + # this test seems unfinished - not asserting anything! convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km) instance = self.new_instance_of_an_optional_code(Rebound, convert_nbody) @@ -240,30 +282,11 @@ def test2(self): instance.particles.copy_values_of_all_attributes_to(stars) stars.savepoint() - if HAS_MATPLOTLIB: - figure = pyplot.figure() - plot = figure.add_subplot(1, 1, 1) - - x_points = earth.get_timeline_of_attribute("x") - y_points = earth.get_timeline_of_attribute("y") - - x_points_in_AU = [t_x[1].value_in(units.AU) for t_x in x_points] - y_points_in_AU = [t_x1[1].value_in(units.AU) for t_x1 in y_points] - - plot.scatter(x_points_in_AU, y_points_in_AU, color="b", marker='o') - - plot.set_xlim(-1.5, 1.5) - plot.set_ylim(-1.5, 1.5) - - test_results_path = self.get_path_to_results() - output_file = os.path.join(test_results_path, "rebound-earth-sun2.svg") - figure.savefig(output_file) - instance.cleanup_code() instance.stop() - def test3(self): - particles = datamodel.Particles(7) + def test_colliding_particles(self): + particles = Particles(7) particles.mass = 0.001 | nbody_system.mass particles.radius = 0.1 | nbody_system.length particles.x = [-101.0, -100.0, -0.5, 0.5, 100.0, 101.0, 104.0] | nbody_system.length @@ -289,7 +312,7 @@ def test3(self): (collisions.particles(0).radius + collisions.particles(1).radius), [True, True, True]) - sticky_merged = datamodel.Particles(len(collisions.particles(0))) + sticky_merged = Particles(len(collisions.particles(0))) sticky_merged.mass = collisions.particles(0).mass + collisions.particles(1).mass sticky_merged.radius = collisions.particles(0).radius for p1, p2, merged in zip(collisions.particles(0), collisions.particles(1), sticky_merged): @@ -315,7 +338,7 @@ def test3(self): [True]) instance.stop() - def test4(self): + def test_energies(self): convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km) instance = self.new_instance_of_an_optional_code(Rebound, convert_nbody) instance.initialize_code() @@ -330,11 +353,11 @@ def test4(self): instance.stop() - def test5(self): + def test_out_of_box(self): instance = self.new_instance_of_an_optional_code(Rebound) instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2 - particles = datamodel.Particles(2) + particles = Particles(2) particles.position = ([0, 0, 0], [1, 0, 0]) | nbody_system.length particles.velocity = ([-1, 0, 0], [2, 0, 0]) | nbody_system.speed particles.radius = 0 | nbody_system.length @@ -350,11 +373,11 @@ def test5(self): self.assertEqual(instance.stopping_conditions.out_of_box_detection.particles(0)[0].key, particles[1].key) instance.stop() - def test6(self): + def test_energies_nbody(self): instance = self.new_instance_of_an_optional_code(Rebound) instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2 - particles = datamodel.Particles(2) + particles = Particles(2) particles.position = ([0, 0, 0], [1, 0, 0]) | nbody_system.length particles.velocity = ([-1, 0, 0], [2, 0, 0]) | nbody_system.speed particles.radius = 0 | nbody_system.length @@ -365,7 +388,7 @@ def test6(self): self.assertAlmostRelativeEquals(instance.potential_energy, particles.potential_energy(G=nbody_system.G)) instance.stop() - def test7(self): + def test_energies_subset(self): instance = self.new_instance_of_an_optional_code(Rebound) instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2 subset0 = 0 @@ -373,13 +396,13 @@ def test7(self): print(subset1) self.assertEqual(subset1, 1) - particles = datamodel.Particles(2) + particles = Particles(2) particles.position = ([0, 0, 0], [1, 0, 0]) | nbody_system.length particles.velocity = ([-1, 0, 0], [2, 0, 0]) | nbody_system.speed particles.radius = 0 | nbody_system.length particles.mass = 0.1 | nbody_system.mass - particles2 = datamodel.Particles(2) + particles2 = Particles(2) particles2.position = ([0, 0, 0], [1, 0, 0]) | nbody_system.length particles2.velocity = ([-1, 0, 0], [2, 0, 0]) | nbody_system.speed particles2.radius = 0 | nbody_system.length @@ -403,13 +426,13 @@ def test8(self): print(subset1) self.assertEqual(subset1, 1) - particles = datamodel.Particles(2) + particles = Particles(2) particles.position = ([0, 0, 0], [1, 0, 0]) | nbody_system.length particles.velocity = ([0, 0, 0], [0, 0.5, 0]) | nbody_system.speed particles.radius = 0 | nbody_system.length particles.mass = 0.1 | nbody_system.mass - particles2 = datamodel.Particles(2) + particles2 = Particles(2) particles2.position = ([0, 0, 0], [2, 0, 0]) | nbody_system.length particles2.velocity = ([0, 0, 0], [0, 1, 0]) | nbody_system.speed particles2.radius = 0 | nbody_system.length @@ -457,8 +480,8 @@ def test9(self): instance = self.new_instance_of_an_optional_code(Rebound) # instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2 - particles = datamodel.Particles(10) - particles.x = numpy.arange(0, 10, 1) | nbody_system.length + particles = Particles(10) + particles.x = np.arange(0, 10, 1) | nbody_system.length particles.y = 0 | nbody_system.length particles.z = 0 | nbody_system.length particles.velocity = [0, 1, 0] | nbody_system.speed @@ -486,7 +509,7 @@ def test10(self): instance = self.new_instance_of_an_optional_code(Rebound) # instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2 - particles = datamodel.Particles(1) + particles = Particles(1) particles.x = 0 | nbody_system.length particles.y = 0 | nbody_system.length particles.z = 0 | nbody_system.length