diff --git a/src/amuse/community/rebound/Makefile b/src/amuse/community/rebound/Makefile index e5e9abcaa9..c8798d3a97 100644 --- a/src/amuse/community/rebound/Makefile +++ b/src/amuse/community/rebound/Makefile @@ -20,15 +20,37 @@ OBJS = interface.o CODELIB = src/librebound.a CODEDIR = src/rebound -LIB_OBJ = $(CODEDIR)/src/rebound.o $(CODEDIR)/src/tree.o $(CODEDIR)/src/particle.o $(CODEDIR)/src/gravity.o \ - $(CODEDIR)/src/integrator.o $(CODEDIR)/src/integrator_whfast.o $(CODEDIR)/src/integrator_whfasthelio.o \ - $(CODEDIR)/src/integrator_ias15.o $(CODEDIR)/src/integrator_sei.o $(CODEDIR)/src/integrator_leapfrog.o \ - $(CODEDIR)/src/integrator_hermes.o $(CODEDIR)/src/boundary.o \ - $(CODEDIR)/src/collision.o $(CODEDIR)/src/tools.o \ - $(CODEDIR)/src/communication_mpi.o $(CODEDIR)/src/display.o $(CODEDIR)/src/derivatives.o \ - $(CODEDIR)/src/glad.o $(CODEDIR)/src/integrator_janus.o $(CODEDIR)/src/transformations.o \ - $(CODEDIR)/src/simulationarchive.o $(CODEDIR)/src/output.o \ - $(CODEDIR)/src/input.o \ +LIB_OBJ = \ + $(CODEDIR)/src/rebound.o \ + $(CODEDIR)/src/tree.o \ + $(CODEDIR)/src/particle.o \ + $(CODEDIR)/src/gravity.o \ + $(CODEDIR)/src/integrator.o \ + $(CODEDIR)/src/integrator_whfast.o \ + $(CODEDIR)/src/integrator_whfast512.o \ + $(CODEDIR)/src/integrator_saba.o \ + $(CODEDIR)/src/integrator_ias15.o \ + $(CODEDIR)/src/integrator_sei.o \ + $(CODEDIR)/src/integrator_bs.o \ + $(CODEDIR)/src/integrator_leapfrog.o \ + $(CODEDIR)/src/integrator_mercurius.o \ + $(CODEDIR)/src/integrator_eos.o \ + $(CODEDIR)/src/boundary.o \ + $(CODEDIR)/src/input.o \ + $(CODEDIR)/src/binarydiff.o \ + $(CODEDIR)/src/output.o \ + $(CODEDIR)/src/collision.o \ + $(CODEDIR)/src/communication_mpi.o \ + $(CODEDIR)/src/display.o \ + $(CODEDIR)/src/tools.o \ + $(CODEDIR)/src/rotations.o \ + $(CODEDIR)/src/derivatives.o \ + $(CODEDIR)/src/simulationarchive.o \ + $(CODEDIR)/src/glad.o \ + $(CODEDIR)/src/integrator_janus.o \ + $(CODEDIR)/src/transformations.o \ + $(CODEDIR)/src/fmemopen.o \ + $(CODEDIR)/src/server.o DOWNLOAD_FROM_WEB = $(PYTHON) ./download.py PATCH_FILES = $(PYTHON) ./patch_files.py @@ -39,25 +61,8 @@ AM_CFLAGS = -I$(AMUSE_DIR)/lib/amuse_mpi all: rebound_worker -DOWNLOAD_CODES?=1 - -ifdef DOWNLOAD_CODES $(CODEDIR)/Makefile: make -C . download -else -$(CODEDIR)/Makefile: - @echo "" - @echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - @echo "" - @echo "DOWNLOAD_CODES is not set. Rebound will not be downloaded and built." - @echo "If you do want Rebound, set DOWNLOAD_CODES to 1." - @echo "bash> export DOWNLOAD_CODES=1" - @echo "csh> setenv DOWNLOAD_CODES 1" - @echo "" - @echo "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" - @echo "" - @make -s --no-print-directory -C . raise_error -endif download: $(RM) -Rf .pc @@ -85,7 +90,7 @@ distclean: $(RM) -Rf .pc $(CODELIB): $(CODEDIR)/Makefile - make -C $(CODEDIR) all CC=$(CC) + make -C $(CODEDIR) librebound CC=$(CC) ar rv $(CODELIB) $(LIB_OBJ) ranlib $(CODELIB) diff --git a/src/amuse/community/rebound/download.py b/src/amuse/community/rebound/download.py index 8a686caa88..953fc1e7ce 100755 --- a/src/amuse/community/rebound/download.py +++ b/src/amuse/community/rebound/download.py @@ -2,18 +2,22 @@ import subprocess import os -import sys -import time +import argparse import urllib.request import urllib.parse import urllib.error -from optparse import OptionParser +from shutil import which -class GetCodeFromHttp(object): - url_template = "https://github.com/hannorein/rebound/archive/{version}.tar.gz" +class GetCodeFromHttp: filename_template = "{version}.tar.gz" - version = "" + name = ["rebound"] + url_template = [ + "https://github.com/hannorein/rebound/archive/{version}.tar.gz" + ] + version = [ + "", + ] def directory(self): return os.path.abspath(os.path.dirname(__file__)) @@ -21,13 +25,13 @@ def directory(self): def src_directory(self): return os.path.join(self.directory(), "src") - def unpack_downloaded_file(self, filename): - print("unpacking", filename) + def unpack_downloaded_file(self, filename, name, version): + print(f"unpacking {filename}") arguments = ["tar", "-xf"] arguments.append(filename) subprocess.call(arguments, cwd=os.path.join(self.src_directory())) subprocess.call( - ["mv", "rebound-{version}".format(version=self.version), "rebound"], + ["mv", f"{name}-{version}", name], cwd=os.path.join(self.src_directory()), ) print("done") @@ -35,42 +39,54 @@ def unpack_downloaded_file(self, filename): def start(self): if os.path.exists("src"): counter = 0 - while os.path.exists("src.{0}".format(counter)): + while os.path.exists(f"src.{counter}"): counter += 1 if counter > 100: print("too many backup directories") break - os.rename("src", "src.{0}".format(counter)) + os.rename("src", f"src.{counter}") os.mkdir("src") - url = self.url_template.format(version=self.version) - filename = self.filename_template.format(version=self.version) - filepath = os.path.join(self.src_directory(), filename) - print("downloading version", self.version, "from", url, "to", filename) - urllib.request.urlretrieve(url, filepath) - print("downloading finished") - self.unpack_downloaded_file(filename) + for i, url_template in enumerate(self.url_template): + url = url_template.format(version=self.version[i]) + filename = self.filename_template.format(version=self.version[i]) + filepath = os.path.join(self.src_directory(), filename) + print(f"downloading version {self.version[i]} from {url} to {filename}") + if which("wget") is not None: + arguments = ["wget", url] + subprocess.call(arguments, cwd=os.path.join(self.src_directory())) + elif which("curl") is not None: + arguments = ["curl", "-L", "-O", url] + subprocess.call(arguments, cwd=os.path.join(self.src_directory())) + else: + urllib.request.urlretrieve(url, filepath) + print("downloading finished") + self.unpack_downloaded_file(filename, self.name[i], self.version[i]) + + +def new_argument_parser(): + result = argparse.ArgumentParser() + result.add_argument( + "--rebound-version", + default="111bcd50de54e4cc70e1d4915918226a4adbe188", + # default="4df2b4603058e7d767aee9c26b586b55a5f4de65", + dest="rebound_version", + help="Rebound commit hash to download", + type=str, + ) + return result def main(version=""): + arguments = new_argument_parser().parse_args() + version = [ + arguments.rebound_version, + ] instance = GetCodeFromHttp() instance.version = version instance.start() -def new_option_parser(): - result = OptionParser() - result.add_option( - "--version", - default="31d117bdc92182073d0941c331f76e95f515bfc6", - dest="version", - help="version number to download", - type="string", - ) - return result - - if __name__ == "__main__": - options, arguments = new_option_parser().parse_args() - main(**options.__dict__) + main() diff --git a/src/amuse/community/rebound/interface.cc b/src/amuse/community/rebound/interface.cc index d9f16247fa..ee9426bd20 100644 --- a/src/amuse/community/rebound/interface.cc +++ b/src/amuse/community/rebound/interface.cc @@ -71,7 +71,7 @@ static inline particle_location get_particle_from_identity(int index_of_the_part for( ReboundSimulationVector::iterator i = codes.begin(); i != codes.end(); i++) { code_state cs = *i; particle.code = cs.code; - particle.p = reb_get_particle_by_hash(particle.code, index_of_the_particle); + particle.p = reb_simulation_particle_by_hash(particle.code, index_of_the_particle); if (particle.p != NULL) break; //*i = cs; } @@ -124,7 +124,7 @@ int set_mass(int index_of_the_particle, double mass){ if(p->m==0){ if(mass>0){ - int index_old=reb_get_particle_index(p); + int index_old=reb_simulation_particle_index(p); if(index_old!=code->N_active){ struct reb_particle tmp = code->particles[index_old]; for(int j=index_old; j>code->N_active; j--){ @@ -137,7 +137,7 @@ int set_mass(int index_of_the_particle, double mass){ } else { if(mass==0){ - int index_old=reb_get_particle_index(p); + int index_old=reb_simulation_particle_index(p); code->N_active--; if(index_old!=code->N_active){ @@ -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; @@ -178,7 +183,7 @@ int new_particle(int * index_of_the_particle, double mass, double x, pt.m = mass; pt.r = radius; pt.hash = new_hash; - reb_add(codes[code_index].code, pt); + reb_simulation_add(codes[code_index].code, pt); //std::cout<<"new particle :"<t+code->dt)*dtsign>=tmax*dtsign && exact_finish_time==1){ - reb_integrator_synchronize(code); + reb_simulation_synchronize(code); code->dt = tmax-code->t; last_step++; }else{ @@ -383,7 +388,7 @@ int _evolve_code(double _tmax, code_state * cs){ } } } - reb_integrator_synchronize(code); + reb_simulation_synchronize(code); code->dt = code->dt_last_done; get_kinetic_energy(cs->subset, &ke1); //printf("Code time: %d , %f -> %f (%f,%f)\n",cs->subset , code->t, tmax, ke1, (ke1-ke)/ke); @@ -437,7 +442,7 @@ int _delete_particle(int index_of_the_particle, int code_index){ if(code_index < 0 || code_index >= (signed) codes.size()){ return -10; } - reb_remove_by_hash(codes[code_index].code, index_of_the_particle, keepSorted); + reb_simulation_remove_particle_by_hash(codes[code_index].code, index_of_the_particle, keepSorted); return 0; } @@ -495,7 +500,7 @@ int get_state(int index_of_the_particle, double * mass, double * x, #endif // COLLISIONS_NONE for( ReboundSimulationVector::iterator i = codes.begin(); i != codes.end(); i++) { code_state cs = *i; - p = reb_get_particle_by_hash(cs.code, index_of_the_particle); + p = reb_simulation_particle_by_hash(cs.code, index_of_the_particle); if (p != NULL) { *subset = cs.subset; break; @@ -608,7 +613,7 @@ int get_subset(int index_of_the_particle, int * subset){ struct reb_particle* p=NULL; for( ReboundSimulationVector::iterator i = codes.begin(); i != codes.end(); i++) { code_state cs = *i; - p = reb_get_particle_by_hash(cs.code, index_of_the_particle); + p = reb_simulation_particle_by_hash(cs.code, index_of_the_particle); if (p != NULL) { *subset = cs.subset; break; @@ -638,7 +643,7 @@ int set_subset(int index_of_the_particle, int subset){ struct reb_particle* p=NULL; for( ReboundSimulationVector::iterator i = codes.begin(); i != codes.end(); i++) { code_state cs = *i; - p = reb_get_particle_by_hash(cs.code, index_of_the_particle); + p = reb_simulation_particle_by_hash(cs.code, index_of_the_particle); if (p != NULL) { if(cs.subset != subset) {return -2;} break; @@ -653,8 +658,8 @@ int cleanup_code() { for( ReboundSimulationVector::iterator i = codes.begin(); i != codes.end(); i++) { code_state cs = *i; if(cs.code){ - reb_remove_all(cs.code); - reb_free_simulation(cs.code); + reb_simulation_remove_all_particles(cs.code); + reb_simulation_free(cs.code); cs.code = 0; *i = cs; } @@ -679,7 +684,7 @@ int initialize_code(){ int nt = omp_get_max_threads(); omp_set_num_threads(nt); #endif - reb_simulation * code = reb_create_simulation(); + reb_simulation * code = reb_simulation_create(); codes.push_back(code_state(code)); code->integrator = reb_simulation::REB_INTEGRATOR_IAS15; code->N_active = 0; @@ -789,8 +794,8 @@ int set_velocity(int index_of_the_particle, double vx, double vy, } int new_subset(int * index, double time_offset) { - reb_simulation * code = reb_create_simulation(); - reb_integrator_reset(code); + reb_simulation * code = reb_simulation_create(); + reb_simulation_reset_integrator(code); code->dt = timestep; if(time_offset < 0) {time_offset = _time;} code->integrator = reb_simulation::REB_INTEGRATOR_IAS15; @@ -811,8 +816,8 @@ int stop_subset(int code_index) { code_state cs = codes[code_index]; if(cs.code) { reb_simulation * code = cs.code; - reb_remove_all(code); - reb_free_simulation(code); + reb_simulation_remove_all_particles(code); + reb_simulation_free(code); cs.code = 0; codes[code_index] = cs; } @@ -845,10 +850,12 @@ int _set_integrator(int value, int code_index){ code->integrator = reb_simulation::REB_INTEGRATOR_LEAPFROG; break; case 5: - code->integrator = reb_simulation::REB_INTEGRATOR_HERMES; + // This integrator was removed + return -1; break; case 6: - code->integrator = reb_simulation::REB_INTEGRATOR_WHFASTHELIO; + // This integrator was removed + return -1; break; case 7: code->integrator = reb_simulation::REB_INTEGRATOR_NONE; @@ -856,6 +863,21 @@ int _set_integrator(int value, int code_index){ case 8: code->integrator = reb_simulation::REB_INTEGRATOR_JANUS; break; + case 9: + 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_EOS; + break; + 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; @@ -1014,7 +1036,7 @@ int set_boundary_size(double boundary_size, int code_index){ return -11; } reb_simulation * code = codes[code_index].code; - reb_configure_box(code,boundary_size,1,1,1); + reb_simulation_configure_box(code,boundary_size,1,1,1); return 0; } diff --git a/src/amuse/community/rebound/interface.py b/src/amuse/community/rebound/interface.py index 0ce390818d..d167c6a796 100644 --- a/src/amuse/community/rebound/interface.py +++ b/src/amuse/community/rebound/interface.py @@ -1,8 +1,17 @@ -from amuse.community import * -from amuse.community.interface.gd import GravitationalDynamicsInterface -from amuse.community.interface.gd import GravitationalDynamics -from amuse.community.interface.gd import SinglePointGravityFieldInterface -from amuse.community.interface.gd import GravityFieldCode +from amuse.units import nbody_system +from amuse.community import ( + CodeInterface, + LiteratureReferencesMixIn, + StoppingConditionInterface, + StoppingConditions, + legacy_function, + LegacyFunctionSpecification, +) +from amuse.community.interface.gd import ( + GravitationalDynamicsInterface, + GravitationalDynamics, + GravityFieldCode, +) class ReboundInterface( @@ -10,7 +19,6 @@ class ReboundInterface( LiteratureReferencesMixIn, GravitationalDynamicsInterface, StoppingConditionInterface, - # SinglePointGravityFieldInterface ): """ REBOUND - An open-source multi-purpose N-body code @@ -105,7 +113,10 @@ def new_particle(): "subset", dtype="int32", direction=function.IN, - description="The subset index of the particle (defaults to 0, use new_subset for higher indices)", + description=( + "The subset index of the particle " + "(defaults to 0, use new_subset for higher indices)" + ), default=0, ) function.result_type = "int32" @@ -180,10 +191,15 @@ def _get_integrator(): "whfast": 1, "sei": 2, "leapfrog": 4, - "hermes": 5, - "whfast-helio": 6, + "hermes": 5, # removed + "whfast-helio": 6, # removed "none": 7, "janus": 8, + "mercurius": 9, + "saba": 10, + "eos": 11, + "bs": 12, + "whfast512": 21, } def set_integrator(self, name, code_index=0): @@ -226,7 +242,12 @@ def _get_solver(): function.can_handle_array = False return function - SOLVERS = {"none": 0, "basic": 1, "compensated": 2, "tree": 3} + SOLVERS = { + "none": 0, + "basic": 1, + "compensated": 2, + "tree": 3, + } def set_solver(self, name, code_index=0): return self._set_solver(self.SOLVERS[name], code_index) @@ -280,7 +301,8 @@ def set_opening_angle2(): def get_eps2(): function = LegacyFunctionSpecification() """ - Get epsilon^2, a softening parameter for gravitational potentials with point particles. + Get epsilon^2, a softening parameter for gravitational potentials with + point particles. """ function = LegacyFunctionSpecification() function.addParameter( @@ -294,7 +316,10 @@ def get_eps2(): "epsilon_squared", dtype="float64", direction=function.OUT, - description="epsilon^2, a softening parameter for gravitational potentials with point particles", + description=( + "epsilon^2, a softening parameter for gravitational " + "potentials with point particles" + ), unit=nbody_system.length * nbody_system.length, ) function.result_type = "int32" @@ -309,14 +334,18 @@ def get_eps2(): @legacy_function def set_eps2(): """ - Set epsilon^2, a softening parameter for gravitational potentials with point particles. + Set epsilon^2, a softening parameter for gravitational potentials with + point particles. """ function = LegacyFunctionSpecification() function.addParameter( "epsilon_squared", dtype="float64", direction=function.IN, - description="epsilon^2, a softening parameter for gravitational potentials with point particles", + description=( + "epsilon^2, a softening parameter for gravitational " + "potentials with point particles" + ), unit=nbody_system.length * nbody_system.length, ) function.addParameter( @@ -338,7 +367,9 @@ def set_eps2(): @legacy_function def _set_boundary(): function = LegacyFunctionSpecification() - function.addParameter("boundary_name", dtype="i", direction=function.IN) + function.addParameter( + "boundary_name", dtype="i", direction=function.IN + ) function.addParameter( "code_index", dtype="int32", @@ -360,12 +391,19 @@ def _get_boundary(): description="Index of the code in rebound", default=0, ) - function.addParameter("boundary_name", dtype="i", direction=function.OUT) + function.addParameter( + "boundary_name", dtype="i", direction=function.OUT + ) function.result_type = "int32" function.can_handle_array = False return function - BOUNDARIES = {"none": 0, "open": 1, "periodic": 2, "shear": 3} + BOUNDARIES = { + "none": 0, + "open": 1, + "periodic": 2, + "shear": 3, + } def set_boundary(self, name, code_index=0): return self._set_boundary(self.BOUNDARIES[name], code_index) @@ -443,7 +481,8 @@ def set_time_step(): """ function = LegacyFunctionSpecification() function.addParameter( - "timestep", dtype="float64", direction=function.IN, description="timestep" + "timestep", dtype="float64", direction=function.IN, + description="timestep" ) function.addParameter( "code_index", @@ -518,7 +557,8 @@ def get_kinetic_energy(): @legacy_function def evolve_model(): """ - Evolve the model until the given time, or until a stopping condition is set. + Evolve the model until the given time, or until a stopping condition is + set. """ function = LegacyFunctionSpecification() function.addParameter( @@ -532,7 +572,9 @@ def evolve_model(): "code_index", dtype="int32", direction=function.IN, - description="Index of the code in rebound (default -1, evolve all systems)", + description=( + "Index of the code in rebound (default -1, evolve all systems)" + ), default=-1, ) function.result_type = "int32" @@ -541,9 +583,9 @@ def evolve_model(): @legacy_function def get_time(): """ - Retrieve the model time. This time should be close to the end time specified - in the evolve code. Or, when a collision was detected, it will be the - model time of the collision. + Retrieve the model time. This time should be close to the end time + specified in the evolve code. Or, when a collision was detected, it + will be the model time of the collision. """ function = LegacyFunctionSpecification() function.addParameter( @@ -599,7 +641,8 @@ def get_time_step(): @legacy_function def new_subset(): """ - Create a new particle subset (and corresponding code). This subset will evolve seperately from others. + Create a new particle subset (and corresponding code). This subset will + evolve separately from others. """ function = LegacyFunctionSpecification() function.can_handle_array = True @@ -616,7 +659,9 @@ def new_subset(): "time_offset", dtype="float64", direction=function.IN, - description="Time of the system (defaults to the current model time)", + description=( + "Time of the system (defaults to the current model time)" + ), default=-1, ) function.result_type = "int32" @@ -629,8 +674,9 @@ def new_subset(): @legacy_function def get_state(): """ - Retrieve the current state of a particle. The *minimal* information of a stellar - dynamics particle (mass, radius, position and velocity) is returned. + Retrieve the current state of a particle. The *minimal* information of + a stellar dynamics particle (mass, radius, position and velocity) is + returned. """ function = LegacyFunctionSpecification() function.can_handle_array = True @@ -638,7 +684,10 @@ def get_state(): "index_of_the_particle", dtype="int32", direction=function.IN, - description="Index of the particle to get the state from. This index must have been returned by an earlier call to :meth:`new_particle`", + description=( + "Index of the particle to get the state from. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), ) function.addParameter( "mass", @@ -713,7 +762,10 @@ def get_subset(): "index_of_the_particle", dtype="int32", direction=function.IN, - description="Index of the particle to get the subset of. This index must have been returned by an earlier call to :meth:`new_particle`", + description=( + "Index of the particle to get the subset of. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), ) function.addParameter( "subset", @@ -725,7 +777,7 @@ def get_subset(): function.can_handle_array = True function.result_doc = """ 0 - OK - particle was found in the model and the information was retreived + particle was found in the model and the information was retrieved -1 - ERROR particle could not be found """ @@ -764,19 +816,25 @@ def set_subset(): "index_of_the_particle", dtype="int32", direction=function.IN, - description="Index of the particle to get the subset of. This index must have been returned by an earlier call to :meth:`new_particle`", + description=( + "Index of the particle to get the subset of. This index must " + "have been returned by an earlier call to :meth:`new_particle`" + ), ) function.addParameter( "subset", dtype="int32", direction=function.IN, - description="The new subset of the particle, as this is actually read only this will fail if changed!", + description=( + "The new subset of the particle, as this is actually read only " + "this will fail if changed!" + ), ) function.result_type = "int32" function.can_handle_array = True function.result_doc = """ 0 - OK - particle was found in the model and the information was retreived + particle was found in the model and the information was retrieved -1 - ERROR particle could not be found """ @@ -828,8 +886,9 @@ def define_parameters(self, handler): "get_solver", "set_solver", "solver", - "name of the gravity solver to use ({0})".format( - sorted(self.SOLVERS.keys()) + ( + f"name of the gravity solver to use " + f"({sorted(self.SOLVERS.keys())})" ), default_value="compensated", ) @@ -846,7 +905,10 @@ def define_parameters(self, handler): "get_opening_angle2", "set_opening_angle2", "opening_angle2", - "opening angle, theta, for building the tree in case of tree solver: between 0 and 1", + ( + "opening angle, theta, for building the tree in case of " + "tree solver: between 0 and 1" + ), default_value=0.5, ) @@ -854,8 +916,9 @@ def define_parameters(self, handler): "get_boundary", "set_boundary", "boundary", - "name of the boundary type to use ({0}) (required for tree solver)".format( - sorted(self.BOUNDARIES.keys()) + ( + f"name of the boundary type to use " + f"({sorted(self.BOUNDARIES.keys())}) (required for tree solver)" ), default_value="none", ) @@ -907,7 +970,9 @@ def define_methods(self, handler): ), ) handler.add_method( - "evolve_model", (nbody_system.time, handler.INDEX), (handler.ERROR_CODE,) + "evolve_model", + (nbody_system.time, handler.INDEX), + (handler.ERROR_CODE,) ) handler.add_method( "get_time", @@ -953,7 +1018,9 @@ def define_methods(self, handler): ), ) handler.add_method( - "get_subset", (handler.NO_UNIT,), (handler.NO_UNIT, handler.ERROR_CODE) + "get_subset", + (handler.NO_UNIT,), + (handler.NO_UNIT, handler.ERROR_CODE), ) handler.add_method( "set_subset", diff --git a/src/amuse/test/suite/codes_tests/test_rebound.py b/src/amuse/test/suite/codes_tests/test_rebound.py index da5933259d..2009f51eb3 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,59 +90,109 @@ 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)) - for result, expected in zip(list(instance.get_position(0).values()), [-0.5, 0.0, 0.0, 0]): + 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]): + 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)) - for result, expected in zip(list(instance.get_position(0).values()), [0.5, 0.0, 0.0, 0]): + 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]): + 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.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 +205,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() @@ -169,28 +227,34 @@ def test8(self): class TestRebound(TestWithMPI): - def new_system_of_sun_and_earth(self): - stars = datamodel.Stars(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.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))) - - return stars - - def test1(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() - interface.parameters.epsilon_squared = 0.0 | units.AU**2 + def evolve_earth_orbit( + self, + integrator=None, + ): + """ + Evolve Earth's orbit, with chosen integrator and other settings + """ + 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(redirection="none") + interface.parameters.epsilon_squared = 0.0 | units.au**2 + if integrator is not None: + interface.parameters.integrator = integrator + self.assertEqual( + interface.parameters.integrator, + integrator, + ) + else: + self.assertEqual( + interface.parameters.integrator, + "ias15", + ) interface.dt_dia = 5000 stars = self.new_system_of_sun_and_earth() @@ -198,37 +262,80 @@ def test1(self): interface.particles.add_particles(stars) - interface.evolve_model(365.0 | units.day) + interface.evolve_model(365.25 | units.day) interface.particles.copy_values_of_all_attributes_to(stars) - position_at_start = earth.position.value_in(units.AU)[0] - position_after_full_rotation = earth.position.value_in(units.AU)[0] - self.assertAlmostRelativeEqual(position_at_start, position_after_full_rotation, 6) + position_at_start = earth.position.value_in(units.au)[0] + position_after_full_rotation = earth.position.value_in(units.au)[0] + self.assertAlmostRelativeEqual( + position_at_start, position_after_full_rotation, 6 + ) - interface.evolve_model(365.0 + (365.0 / 2) | units.day) + interface.evolve_model(365.25 + (365.25 / 2) | units.day) - self.assertAlmostRelativeEqual(interface.model_time, 365.0 + (365.0 / 2) | units.day, 3) + self.assertAlmostRelativeEqual( + interface.model_time, + 365.25 + (365.25 / 2) | units.day, 3 + ) interface.particles.copy_values_of_all_attributes_to(stars) - position_after_half_a_rotation = earth.position.value_in(units.AU)[0] - self.assertAlmostRelativeEqual(-position_at_start, position_after_half_a_rotation, 3) - - interface.evolve_model(365.0 + (365.0 / 2) + (365.0 / 4) | units.day) - self.assertAlmostRelativeEqual(interface.model_time, 365.0 + (365.0 / 2) + (365.0 / 4) | units.day, 3) + position_after_half_a_rotation = earth.position.value_in(units.au)[0] + self.assertAlmostRelativeEqual( + -position_at_start, position_after_half_a_rotation, 3 + ) + + interface.evolve_model( + 365.25 + (365.25 / 2) + (365.25 / 4) | units.day + ) + self.assertAlmostRelativeEqual( + interface.model_time, + 365.25 + (365.25 / 2) + (365.25 / 4) | units.day, 3 + ) interface.particles.copy_values_of_all_attributes_to(stars) - position_after_half_a_rotation = earth.position.value_in(units.AU)[1] - self.assertAlmostRelativeEqual(-position_at_start, position_after_half_a_rotation, 3) + position_after_half_a_rotation = earth.position.value_in(units.au)[1] + self.assertAlmostRelativeEqual( + -position_at_start, position_after_half_a_rotation, 3 + ) interface.cleanup_code() - interface.stop() - def test2(self): - convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km) + def new_system_of_sun_and_earth(self): + stars = Particles(2) + sun = stars[0] + sun.mass = units.MSun(1.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(np.array((149.5e6, 0.0, 0.0))) + earth.velocity = units.ms(np.array((0.0, 29800, 0.0))) + + return stars + + def test_evolve_default_integrator(self): + self.evolve_earth_orbit() - instance = self.new_instance_of_an_optional_code(Rebound, convert_nbody) + def test_evolve_ias15_integrator(self): + self.evolve_earth_orbit(integrator="ias15") + + def test_evolve_whfast_integrator(self): + self.evolve_earth_orbit(integrator="whfast") + + def xtest_plot(self): + # this test was meant to plot - but that's not a proper test + 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() - instance.parameters.epsilon_squared = 0.0 | units.AU**2 + instance.parameters.epsilon_squared = 0.0 | units.au**2 instance.dt_dia = 5000 stars = self.new_system_of_sun_and_earth() @@ -240,36 +347,21 @@ 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 + particles.x = [ + -101.0, -100.0, -0.5, 0.5, 100.0, 101.0, 104.0 + ] | nbody_system.length particles.y = 0 | nbody_system.length particles.z = 0 | nbody_system.length - particles.velocity = [[2, 0, 0], [-2, 0, 0]]*3 + [[-4, 0, 0]] | nbody_system.speed + particles.velocity = ( + [[2, 0, 0], [-2, 0, 0]]*3 + [[-4, 0, 0]] + ) | nbody_system.speed instance = self.new_instance_of_an_optional_code(Rebound) instance.initialize_code() @@ -284,21 +376,40 @@ def test3(self): self.assertTrue(instance.model_time < 0.5 | nbody_system.time) self.assertEqual(len(collisions.particles(0)), 3) self.assertEqual(len(collisions.particles(1)), 3) - self.assertEqual(len(particles - collisions.particles(0) - collisions.particles(1)), 1) - self.assertEqual(abs(collisions.particles(0).x - collisions.particles(1).x) <= - (collisions.particles(0).radius + collisions.particles(1).radius), - [True, True, True]) - - sticky_merged = datamodel.Particles(len(collisions.particles(0))) - sticky_merged.mass = collisions.particles(0).mass + collisions.particles(1).mass + self.assertEqual( + len( + particles - collisions.particles(0) - collisions.particles(1) + ), + 1 + ) + self.assertEqual( + abs(collisions.particles(0).x - collisions.particles(1).x) + <= ( + collisions.particles(0).radius + + collisions.particles(1).radius + ), + [True, True, True] + ) + + 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): + for p1, p2, merged in zip( + collisions.particles(0), + collisions.particles(1), + sticky_merged + ): merged.position = (p1 + p2).center_of_mass() merged.velocity = (p1 + p2).center_of_mass_velocity() print(instance.model_time) print(instance.particles) - instance.particles.remove_particles(collisions.particles(0) + collisions.particles(1)) + instance.particles.remove_particles( + collisions.particles(0) + + collisions.particles(1) + ) instance.particles.add_particles(sticky_merged) print(instance.particles) instance.evolve_model(1.0 | nbody_system.time) @@ -309,32 +420,51 @@ def test3(self): self.assertTrue(instance.model_time < 1.0 | nbody_system.time) self.assertEqual(len(collisions.particles(0)), 1) self.assertEqual(len(collisions.particles(1)), 1) - self.assertEqual(len(instance.particles - collisions.particles(0) - collisions.particles(1)), 2) - self.assertEqual(abs(collisions.particles(0).x - collisions.particles(1).x) <= - (collisions.particles(0).radius + collisions.particles(1).radius), - [True]) + self.assertEqual( + len( + instance.particles + - collisions.particles(0) + - collisions.particles(1) + ), 2 + ) + self.assertEqual( + abs(collisions.particles(0).x - collisions.particles(1).x) + <= ( + collisions.particles(0).radius + + collisions.particles(1).radius + ), + [True] + ) instance.stop() - def test4(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) + 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() - instance.parameters.epsilon_squared = 0.0 | units.AU**2 + instance.parameters.epsilon_squared = 0.0 | units.au**2 instance.dt_dia = 5000 stars = self.new_system_of_sun_and_earth() earth = stars[1] instance.particles.add_particles(stars) - self.assertAlmostRelativeEquals(instance.kinetic_energy, stars.kinetic_energy()) - self.assertAlmostRelativeEquals(instance.potential_energy, stars.potential_energy(), 10) + self.assertAlmostRelativeEquals( + instance.kinetic_energy, stars.kinetic_energy() + ) + self.assertAlmostRelativeEquals( + instance.potential_energy, stars.potential_energy(), 10 + ) 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 @@ -345,27 +475,41 @@ def test5(self): instance.parameters.stopping_conditions_out_of_box_size = 2 | nbody_system.length instance.parameters.stopping_conditions_out_of_box_use_center_of_mass = False instance.evolve_model(1 | nbody_system.time) - self.assertTrue(instance.stopping_conditions.out_of_box_detection.is_set()) - self.assertEqual(len(instance.stopping_conditions.out_of_box_detection.particles(0)), 1) - self.assertEqual(instance.stopping_conditions.out_of_box_detection.particles(0)[0].key, particles[1].key) + self.assertTrue( + instance.stopping_conditions.out_of_box_detection.is_set() + ) + self.assertEqual( + len( + instance.stopping_conditions.out_of_box_detection.particles(0) + ), 1 + ) + 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 particles.mass = 0.1 | nbody_system.mass instance.particles.add_particles(particles) - self.assertAlmostRelativeEquals(instance.kinetic_energy, particles.kinetic_energy()) - self.assertAlmostRelativeEquals(instance.potential_energy, particles.potential_energy(G=nbody_system.G)) + self.assertAlmostRelativeEquals( + instance.kinetic_energy, particles.kinetic_energy() + ) + 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 +517,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 @@ -389,10 +533,21 @@ def test7(self): instance.particles.add_particles(particles) instance.particles.add_particles(particles2) print(instance.particles) - self.assertAlmostRelativeEquals(instance.kinetic_energy, particles.kinetic_energy()) - self.assertAlmostRelativeEquals(instance.potential_energy, particles.potential_energy(G=nbody_system.G)) - self.assertAlmostRelativeEquals(instance.get_kinetic_energy(subset1), particles2.kinetic_energy()) - self.assertAlmostRelativeEquals(instance.get_potential_energy(subset1), particles2.potential_energy(G=nbody_system.G)) + self.assertAlmostRelativeEquals( + instance.kinetic_energy, particles.kinetic_energy() + ) + self.assertAlmostRelativeEquals( + instance.potential_energy, + particles.potential_energy(G=nbody_system.G) + ) + self.assertAlmostRelativeEquals( + instance.get_kinetic_energy(subset1), + particles2.kinetic_energy() + ) + self.assertAlmostRelativeEquals( + instance.get_potential_energy(subset1), + particles2.potential_energy(G=nbody_system.G) + ) # instance.stop() def test8(self): @@ -403,13 +558,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 @@ -419,10 +574,24 @@ def test8(self): instance.particles.add_particles(particles) instance.particles.add_particles(particles2) instance.evolve_model(0.1 | nbody_system.time) - self.assertAlmostRelativeEquals(instance.kinetic_energy, particles.kinetic_energy(), 2) - self.assertAlmostRelativeEquals(instance.potential_energy, particles.potential_energy(G=nbody_system.G), 2) - self.assertAlmostRelativeEquals(instance.get_kinetic_energy(subset1), particles2.kinetic_energy(), 2) - self.assertAlmostRelativeEquals(instance.get_potential_energy(subset1), particles2.potential_energy(G=nbody_system.G), 2) + self.assertAlmostRelativeEquals( + instance.kinetic_energy, particles.kinetic_energy(), 2 + ) + self.assertAlmostRelativeEquals( + instance.potential_energy, + particles.potential_energy(G=nbody_system.G), + 2 + ) + self.assertAlmostRelativeEquals( + instance.get_kinetic_energy(subset1), + particles2.kinetic_energy(), + 2 + ) + self.assertAlmostRelativeEquals( + instance.get_potential_energy(subset1), + particles2.potential_energy(G=nbody_system.G), + 2 + ) particles_evolved = instance.particles.copy() instance.stop() @@ -435,8 +604,16 @@ def test8(self): print("HHH:", particles_evolved.subset) print("p2e:", particles1evolved) print(instance1.particles) - self.assertAlmostRelativeEquals(instance1.particles.position, particles1evolved.position, 10) - self.assertAlmostRelativeEquals(instance1.particles.velocity, particles1evolved.velocity, 10) + self.assertAlmostRelativeEquals( + instance1.particles.position, + particles1evolved.position, + 10 + ) + self.assertAlmostRelativeEquals( + instance1.particles.velocity, + particles1evolved.velocity, + 10 + ) instance1.stop() instance2 = self.new_instance_of_an_optional_code(Rebound) @@ -445,20 +622,30 @@ def test8(self): particles2c.subset = 0 instance2.particles.add_particles(particles2c) instance2.evolve_model(0.1 | nbody_system.time) - particles2evolved = particles_evolved[particles_evolved.subset == subset1] + particles2evolved = particles_evolved[ + particles_evolved.subset == subset1 + ] print("HHH:", particles_evolved.subset) print("p2e:", particles2evolved) print(instance2.particles) - self.assertAlmostRelativeEquals(instance2.particles.position, particles2evolved.position, 10) - self.assertAlmostRelativeEquals(instance2.particles.velocity, particles2evolved.velocity, 10) + self.assertAlmostRelativeEquals( + instance2.particles.position, + particles2evolved.position, + 10 + ) + self.assertAlmostRelativeEquals( + instance2.particles.velocity, + particles2evolved.velocity, + 10 + ) instance2.stop() 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 @@ -474,7 +661,9 @@ def test9(self): self.assertEqual(instance.particles[5].mass, 1 | nbody_system.mass) self.assertEqual(instance.particles[6].mass, 0 | nbody_system.mass) self.assertEqual(instance.particles[5].x, 5 | nbody_system.length) - self.assertAlmostRelativeEquals(instance.particles[6].x, 5.5403023 | nbody_system.length, 5) + self.assertAlmostRelativeEquals( + instance.particles[6].x, 5.5403023 | nbody_system.length, 5 + ) instance.particles.remove_particle(particles[4]) self.assertEqual(instance.particles[4].mass, 1 | nbody_system.mass) self.assertEqual(instance.particles[5].mass, 0 | nbody_system.mass) @@ -486,7 +675,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 @@ -496,6 +685,13 @@ def test10(self): instance.particles.add_particles(particles) instance.evolve_model(4 | nbody_system.time) - self.assertAlmostRelativeEquals(instance.particles[0].velocity, [0, 2, 0] | nbody_system.speed) - self.assertAlmostRelativeEquals(instance.particles[0].position, [0, 2 * 4, 0] | nbody_system.length, 8) + self.assertAlmostRelativeEquals( + instance.particles[0].velocity, + [0, 2, 0] | nbody_system.speed + ) + self.assertAlmostRelativeEquals( + instance.particles[0].position, + [0, 2 * 4, 0] | nbody_system.length, + 8 + ) instance.stop()