From 2211a20f1b8d772544796b6c4e1467410f5be883 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Aug 2024 09:26:06 +0200 Subject: [PATCH 1/6] Standardise class names, keep compatibility --- src/amuse/community/sse/interface.py | 610 ++++++++++++++++----------- 1 file changed, 358 insertions(+), 252 deletions(-) diff --git a/src/amuse/community/sse/interface.py b/src/amuse/community/sse/interface.py index 22fd47b130..0fb35b18e5 100644 --- a/src/amuse/community/sse/interface.py +++ b/src/amuse/community/sse/interface.py @@ -1,157 +1,228 @@ -import numpy +""" +AMUSE interface for the stellar population synthesis code SSE +""" + from operator import itemgetter -from amuse.community import * -from amuse.units import units -from amuse.units import constants +import numpy +from amuse.units import units, constants from amuse.support.interface import InCodeComponentImplementation -from amuse.community.interface import common from amuse.datamodel import Particles from amuse.datamodel import ParticlesSubset -class SSEInterface(CodeInterface, common.CommonCodeInterface , LiteratureReferencesMixIn): + +from amuse.community.interface import common +from amuse.community import ( + CodeInterface, + LiteratureReferencesMixIn, + legacy_function, + LegacyFunctionSpecification, + NO_UNIT, +) + + +class SseInterface( + CodeInterface, common.CommonCodeInterface, LiteratureReferencesMixIn +): """ Stellar evolution is performed by the rapid single-star evolution (SSE) - algorithm. This is a package of analytical formulae fitted to the detailed - models of Pols et al. (1998) that covers all phases of evolution from the - zero-age main-sequence up to and including remnant phases. It is valid for - masses in the range 0.1-100 Msun and metallicity can be varied. The SSE - package contains a prescription for mass loss by stellar winds. It also - follows the evolution of rotational angular momentum for the star. Full + algorithm. This is a package of analytical formulae fitted to the detailed + models of Pols et al. (1998) that covers all phases of evolution from the + zero-age main-sequence up to and including remnant phases. It is valid for + masses in the range 0.1-100 Msun and metallicity can be varied. The SSE + package contains a prescription for mass loss by stellar winds. It also + follows the evolution of rotational angular momentum for the star. Full details can be found in the SSE paper: - + .. [#] ADS:2000MNRAS.315..543H (Hurley J.R., Pols O.R., Tout C.A., 2000, MNRAS, 315, 543: .. [#] ... "Comprehensive analytic formulae for stellar evolution as a function of mass and metallicity") """ + def __init__(self, **options): CodeInterface.__init__(self, name_of_the_worker="sse_worker", **options) LiteratureReferencesMixIn.__init__(self) - - @legacy_function + + @legacy_function def initialize(): - function = LegacyFunctionSpecification() - function.addParameter('z_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('neta_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('bwind_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('hewind_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('sigma_in', dtype='d', direction=function.IN, unit = units.km / units.s) - function.addParameter('ifflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('wdflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('bhflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('nsflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('mxns_in', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('pts1_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('pts2_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('pts3_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('status', dtype='i', direction=function.OUT, unit = NO_UNIT) + function = LegacyFunctionSpecification() + function.addParameter("z_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("neta_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter( + "bwind_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "hewind_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "sigma_in", dtype="d", direction=function.IN, unit=units.km / units.s + ) + function.addParameter( + "ifflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "wdflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "bhflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "nsflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "mxns_in", dtype="d", direction=function.IN, unit=units.MSun + ) + function.addParameter("pts1_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("pts2_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("pts3_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("status", dtype="i", direction=function.OUT, unit=NO_UNIT) return function - - @legacy_function + + @legacy_function def evolve_star(): - function = LegacyFunctionSpecification() - function.name = 'evolve0' - function.can_handle_array = True - function.addParameter('kw', dtype='i', direction=function.INOUT, unit = units.stellar_type) - function.addParameter('mass', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('mt', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('r', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('lum', dtype='d', direction=function.INOUT, unit = units.LSun) - function.addParameter('mc', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('rc', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('menv', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('renv', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('ospin', dtype='d', direction=function.INOUT, unit = units.yr**-1) - function.addParameter('epoch', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('tm', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('tphys', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('tphysf', dtype='d', direction=function.INOUT, unit = units.Myr) + function = LegacyFunctionSpecification() + function.name = "evolve0" + function.can_handle_array = True + function.addParameter( + "kw", dtype="i", direction=function.INOUT, unit=units.stellar_type + ) + function.addParameter( + "mass", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "mt", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter("r", dtype="d", direction=function.INOUT, unit=units.RSun) + function.addParameter( + "lum", dtype="d", direction=function.INOUT, unit=units.LSun + ) + function.addParameter( + "mc", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "rc", dtype="d", direction=function.INOUT, unit=units.RSun + ) + function.addParameter( + "menv", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "renv", dtype="d", direction=function.INOUT, unit=units.RSun + ) + function.addParameter( + "ospin", dtype="d", direction=function.INOUT, unit=units.yr**-1 + ) + function.addParameter( + "epoch", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter("tm", dtype="d", direction=function.INOUT, unit=units.Myr) + function.addParameter( + "tphys", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter( + "tphysf", dtype="d", direction=function.INOUT, unit=units.Myr + ) return function - - @legacy_function + + @legacy_function def get_time_step(): - function = LegacyFunctionSpecification() - function.can_handle_array = True - function.addParameter('kw', dtype='i', direction=function.IN, unit = units.stellar_type) - function.addParameter('mass', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('age', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('mt', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('tm', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('epoch', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('dt', dtype='d', direction=function.OUT, unit = units.Myr) - + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "kw", dtype="i", direction=function.IN, unit=units.stellar_type + ) + function.addParameter("mass", dtype="d", direction=function.IN, unit=units.MSun) + function.addParameter("age", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter("mt", dtype="d", direction=function.IN, unit=units.MSun) + function.addParameter("tm", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter("epoch", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter("dt", dtype="d", direction=function.OUT, unit=units.Myr) + return function - - @legacy_function + + @legacy_function def get_mass_loss_wind(): - function = LegacyFunctionSpecification() - function.can_handle_array = True - function.addParameter('kw', dtype='i', direction=function.IN, unit = units.stellar_type) - function.addParameter('lum', dtype='d', direction=function.IN, unit = units.LSun) - function.addParameter('r', dtype='d', direction=function.IN, unit = units.RSun) - function.addParameter('mt', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('mc', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('mlout', dtype='d', direction=function.OUT, unit = units.MSun/units.yr) - + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "kw", dtype="i", direction=function.IN, unit=units.stellar_type + ) + function.addParameter("lum", dtype="d", direction=function.IN, unit=units.LSun) + function.addParameter("r", dtype="d", direction=function.IN, unit=units.RSun) + function.addParameter("mt", dtype="d", direction=function.IN, unit=units.MSun) + function.addParameter("mc", dtype="d", direction=function.IN, unit=units.MSun) + function.addParameter( + "mlout", dtype="d", direction=function.OUT, unit=units.MSun / units.yr + ) + return function - - @legacy_function + @legacy_function def get_gyration_radius(): - function = LegacyFunctionSpecification() - function.can_handle_array = True - function.addParameter('kw', dtype='i', direction=function.IN, unit = units.stellar_type) - function.addParameter('mass', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('mt', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('r', dtype='d', direction=function.IN, unit = units.RSun) - function.addParameter('lum', dtype='d', direction=function.IN, unit = units.LSun) - function.addParameter('epoch', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('tm', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('tphys', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('rg', dtype='d', direction=function.OUT, unit = NO_UNIT) - + function = LegacyFunctionSpecification() + function.can_handle_array = True + function.addParameter( + "kw", dtype="i", direction=function.IN, unit=units.stellar_type + ) + function.addParameter("mass", dtype="d", direction=function.IN, unit=units.MSun) + function.addParameter("mt", dtype="d", direction=function.IN, unit=units.MSun) + function.addParameter("r", dtype="d", direction=function.IN, unit=units.RSun) + function.addParameter("lum", dtype="d", direction=function.IN, unit=units.LSun) + function.addParameter("epoch", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter("tm", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter("tphys", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter("rg", dtype="d", direction=function.OUT, unit=NO_UNIT) + return function - + def initialize_code(self): return 0 - + def commit_parameters(self): return 0 - + def recommit_parameters(self): return 0 - + def cleanup_code(self): return 0 - + def commit_particles(self): return 0 - - - -class SSEParticles(Particles): - - def __init__(self, code_interface, storage = None): - Particles.__init__(self, storage = storage) - self._private.code_interface = code_interface - self.add_calculated_attribute("temperature", self.calculate_effective_temperature, ["luminosity", "radius"]) - self.add_function_attribute("evolve_one_step", self.particleset_evolve_one_step, self.evolve_one_step) - self.add_function_attribute("evolve_for", self.particleset_evolve_for, self.evolve_for) - + + +class SseParticles(Particles): + + def __init__(self, code_interface, storage=None): + Particles.__init__(self, storage=storage) + self._private.code_interface = code_interface + self.add_calculated_attribute( + "temperature", + self.calculate_effective_temperature, + ["luminosity", "radius"], + ) + self.add_function_attribute( + "evolve_one_step", self.particleset_evolve_one_step, self.evolve_one_step + ) + self.add_function_attribute( + "evolve_for", self.particleset_evolve_for, self.evolve_for + ) + def calculate_effective_temperature(self, luminosity, radius): - return ((luminosity/(constants.four_pi_stefan_boltzmann*radius**2))**.25).in_(units.K) - - def add_particles_to_store(self, keys, attributes = [], values = []): + return ( + (luminosity / (constants.four_pi_stefan_boltzmann * radius**2)) ** 0.25 + ).in_(units.K) + + def add_particles_to_store(self, keys, attributes=[], values=[]): if len(keys) == 0: return - + all_attributes = [] all_attributes.extend(attributes) all_values = [] all_values.extend(values) - + mapping_from_attribute_to_default_value = { - "stellar_type" : 1 | units.stellar_type, + "stellar_type": 1 | units.stellar_type, "radius": 0 | units.RSun, - "luminosity": 0 | units.LSun, + "luminosity": 0 | units.LSun, "core_mass": 0 | units.MSun, "CO_core_mass": 0 | units.MSun, "core_radius": 0 | units.RSun, @@ -160,221 +231,232 @@ def add_particles_to_store(self, keys, attributes = [], values = []): "epoch": 0 | units.Myr, "spin": 0 | units.yr**-1, "main_sequence_lifetime": 0 | units.Myr, - "age": 0 | units.Myr + "age": 0 | units.Myr, } - + given_attributes = set(attributes) - - if not "initial_mass" in given_attributes: + + if "initial_mass" not in given_attributes: index_of_mass_attibute = attributes.index("mass") all_attributes.append("initial_mass") all_values.append(values[index_of_mass_attibute] * 1.0) - + for attribute, default_value in mapping_from_attribute_to_default_value.items(): - if not attribute in given_attributes: + if attribute not in given_attributes: all_attributes.append(attribute) all_values.append(default_value.as_vector_with_length(len(keys))) - - super(SSEParticles, self).add_particles_to_store(keys, all_attributes, all_values) - + + super(SseParticles, self).add_particles_to_store( + keys, all_attributes, all_values + ) + added_particles = ParticlesSubset(self, keys) self._private.code_interface._evolve_particles(added_particles, 0 | units.yr) - + def evolve_one_step(self, particles, subset): - self._private.code_interface._evolve_particles(subset.as_set(), subset.age + subset.time_step) - + self._private.code_interface._evolve_particles( + subset.as_set(), subset.age + subset.time_step + ) + def particleset_evolve_one_step(self, particles): - self._private.code_interface._evolve_particles(particles, particles.age + particles.time_step) - + self._private.code_interface._evolve_particles( + particles, particles.age + particles.time_step + ) + def evolve_for(self, particles, subset, delta_time): - self._private.code_interface._evolve_particles(subset.as_set(), subset.age + delta_time) - + self._private.code_interface._evolve_particles( + subset.as_set(), subset.age + delta_time + ) + def particleset_evolve_for(self, particles, delta_time): - self._private.code_interface._evolve_particles(particles, particles.age + delta_time) - - + self._private.code_interface._evolve_particles( + particles, particles.age + delta_time + ) + def get_defined_attribute_names(self): return ["mass", "radius"] - - - - - - - +class Sse(common.CommonCode): -class SSE(common.CommonCode): - def __init__(self, **options): - InCodeComponentImplementation.__init__(self, SSEInterface(**options), **options) - + InCodeComponentImplementation.__init__(self, SseInterface(**options), **options) + self.model_time = 0.0 | units.yr - - + def define_parameters(self, handler): - + handler.add_caching_parameter( - "initialize", - "z_in", - "metallicity", - "Metallicity of all stars", - 0.02 + "initialize", "z_in", "metallicity", "Metallicity of all stars", 0.02 ) - + handler.add_caching_parameter( "initialize", "neta_in", "reimers_mass_loss_coefficient", "Reimers mass-loss coefficient (neta*4x10^-13; 0.5 normally)", - 0.5 + 0.5, ) - + handler.add_caching_parameter( "initialize", "bwind_in", "binary_enhanced_mass_loss_parameter", "The binary enhanced mass loss parameter (inactive for single).", - 0.0 + 0.0, ) - + handler.add_caching_parameter( "initialize", "hewind_in", "helium_star_mass_loss_factor", "Helium star mass loss factor", - 0.5 + 0.5, ) - + handler.add_caching_parameter( "initialize", "sigma_in", "SN_kick_speed_dispersion", "The dispersion in the Maxwellian for the SN kick speed (190 km/s).", - 190.0 | units.km / units.s + 190.0 | units.km / units.s, ) - + handler.add_caching_parameter( "initialize", "ifflag_in", "white_dwarf_IFMR_flag", - "ifflag > 0 uses white dwarf IFMR (initial-final mass relation) of HPE, 1995, MNRAS, 272, 800 (0).", - 0 + "ifflag > 0 uses white dwarf IFMR (initial-final mass relation) of HPE, " + "1995, MNRAS, 272, 800 (0).", + 0, ) - + handler.add_caching_parameter( "initialize", "wdflag_in", "white_dwarf_cooling_flag", "wdflag > 0 uses modified-Mestel cooling for WDs (0).", - 1 + 1, ) - + handler.add_caching_parameter( "initialize", "bhflag_in", "black_hole_kick_flag", "bhflag > 0 allows velocity kick at BH formation (0).", - 0 + 0, ) - + handler.add_caching_parameter( "initialize", "nsflag_in", "neutron_star_mass_flag", - "nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1).", - 1 + "nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, " + "407 (1).", + 1, ) - + handler.add_caching_parameter( "initialize", "mxns_in", "maximum_neutron_star_mass", "The maximum neutron star mass (1.8, nsflag=0; 3.0, nsflag=1).", - 3.0 | units.MSun + 3.0 | units.MSun, ) - + handler.add_caching_parameter( "initialize", "pts1_in", "fractional_time_step_1", - "The timesteps chosen in each evolution phase as decimal fractions of the time taken in that phase: MS (0.05)", - 0.05 + "The timesteps chosen in each evolution phase as decimal fractions of the " + "time taken in that phase: MS (0.05)", + 0.05, ) - + handler.add_caching_parameter( "initialize", "pts2_in", "fractional_time_step_2", - "The timesteps chosen in each evolution phase as decimal fractions of the time taken in that phase: GB, CHeB, AGB, HeGB (0.01)", - 0.01 + "The timesteps chosen in each evolution phase as decimal fractions of the " + "time taken in that phase: GB, CHeB, AGB, HeGB (0.01)", + 0.01, ) - + handler.add_caching_parameter( "initialize", "pts3_in", "fractional_time_step_3", - "The timesteps chosen in each evolution phase as decimal fractions of the time taken in that phase: HG, HeMS (0.02)", - 0.02 + "The timesteps chosen in each evolution phase as decimal fractions of the " + "time taken in that phase: HG, HeMS (0.02)", + 0.02, ) - + def define_state(self, handler): common.CommonCode.define_state(self, handler) - handler.add_transition('INITIALIZED','RUN','commit_parameters') - handler.add_method('RUN', 'evolve_star') - - handler.add_method('RUN', 'before_get_parameter') - handler.add_method('RUN', 'before_set_parameter') - - + handler.add_transition("INITIALIZED", "RUN", "commit_parameters") + handler.add_method("RUN", "evolve_star") + + handler.add_method("RUN", "before_get_parameter") + handler.add_method("RUN", "before_set_parameter") + def define_particle_sets(self, handler): - handler.define_inmemory_set('particles', SSEParticles) - + handler.define_inmemory_set("particles", SseParticles) + handler.add_attribute( - 'particles', - 'time_step', - 'get_time_step', - ('stellar_type', 'initial_mass', 'age', - 'mass', 'main_sequence_lifetime', 'epoch') + "particles", + "time_step", + "get_time_step", + ( + "stellar_type", + "initial_mass", + "age", + "mass", + "main_sequence_lifetime", + "epoch", + ), ) - + handler.add_attribute( - 'particles', - 'mass_loss_wind', - 'get_mass_loss_wind', - ('stellar_type', 'luminosity', - 'radius', 'mass', - 'CO_core_mass') - ) - + "particles", + "mass_loss_wind", + "get_mass_loss_wind", + ("stellar_type", "luminosity", "radius", "mass", "CO_core_mass"), + ) + handler.add_attribute( - 'particles', - 'gyration_radius', - 'get_gyration_radius', - ('stellar_type', 'initial_mass','mass','radius', - 'luminosity','epoch','main_sequence_lifetime', - 'age') - ) - + "particles", + "gyration_radius", + "get_gyration_radius", + ( + "stellar_type", + "initial_mass", + "mass", + "radius", + "luminosity", + "epoch", + "main_sequence_lifetime", + "age", + ), + ) + def _evolve_particles(self, particles, end_time): attributes = ( "stellar_type", - "initial_mass", - "mass", - "radius", - "luminosity", - "core_mass", + "initial_mass", + "mass", + "radius", + "luminosity", + "core_mass", "core_radius", - "convective_envelope_mass", - "convective_envelope_radius", - "spin", + "convective_envelope_mass", + "convective_envelope_radius", + "spin", "epoch", "main_sequence_lifetime", - "age", - "end_time" + "age", + "end_time", ) - + result = self.evolve_star( particles.stellar_type, particles.initial_mass, @@ -389,65 +471,88 @@ def _evolve_particles(self, particles, end_time): particles.epoch, particles.main_sequence_lifetime, particles.age, - end_time.as_vector_with_length(len(particles))) - + end_time.as_vector_with_length(len(particles)), + ) + # For helium (and helium exhausted) stars, the core mass returned is actually the CO core mass - type = result[0].value_in(units.stellar_type) - helium_star_selection = (type > 6) & (type < 16) & (type != 10) + star_type = result[0].value_in(units.stellar_type) + helium_star_selection = (star_type > 6) & (star_type < 16) & (star_type != 10) helium_stars = particles[helium_star_selection] other_stars = particles - helium_stars if len(helium_stars): helium_stars_results = [sub[helium_star_selection] for sub in result] helium_stars_results.append(helium_stars_results[2]) - helium_stars.set_values_in_store(helium_stars.get_all_indices_in_store(), ( - "stellar_type", "initial_mass", "mass", "radius", "luminosity", - "CO_core_mass", - "core_radius", "convective_envelope_mass", "convective_envelope_radius", "spin", "epoch", - "main_sequence_lifetime", "age", "end_time", - "core_mass"), helium_stars_results) + helium_stars.set_values_in_store( + helium_stars.get_all_indices_in_store(), + ( + "stellar_type", + "initial_mass", + "mass", + "radius", + "luminosity", + "CO_core_mass", + "core_radius", + "convective_envelope_mass", + "convective_envelope_radius", + "spin", + "epoch", + "main_sequence_lifetime", + "age", + "end_time", + "core_mass", + ), + helium_stars_results, + ) if len(other_stars): other_star_selection = numpy.logical_not(helium_star_selection) - other_stars.set_values_in_store(other_stars.get_all_indices_in_store(), attributes, - [sub[other_star_selection] for sub in result]) - - def evolve_model(self, end_time = None, keep_synchronous = True): + other_stars.set_values_in_store( + other_stars.get_all_indices_in_store(), + attributes, + [sub[other_star_selection] for sub in result], + ) + + def evolve_model(self, end_time=None, keep_synchronous=True): if not keep_synchronous: - self._evolve_particles(self.particles, self.particles.time_step + self.particles.age) + self._evolve_particles( + self.particles, self.particles.time_step + self.particles.age + ) return - + if end_time is None: end_time = self.model_time + min(self.particles.time_step) - self._evolve_particles(self.particles, end_time - self.model_time + self.particles.age) + self._evolve_particles( + self.particles, end_time - self.model_time + self.particles.age + ) self.model_time = end_time - - def _evolve_model_old(self, end_time = None, keep_synchronous = True): + + def _evolve_model_old(self, end_time=None, keep_synchronous=True): """ - This is the old implementation of evolve_model. Even with (keep_synchronous = True) - it is unable to evolve all stars to a common age, since it relies on the - individual timesteps as determined by the community code. Furthermore, it - is not suited to simulations with ongoing star formation, since it evolves - newly created stars to the same age as the old stars. + This is the old implementation of evolve_model. Even with + (keep_synchronous = True) it is unable to evolve all stars to a common + age, since it relies on the individual timesteps as determined by the + community code. Furthermore, it is not suited to simulations with + ongoing star formation, since it evolves newly created stars to the + same age as the old stars. """ if end_time is None: if keep_synchronous: ages = self.particles.age index, min_age = min(enumerate(ages), key=itemgetter(1)) new_age = min_age + self.particles[index].time_step - selection = self.particles.select(lambda x : x < new_age, ["age"]) + selection = self.particles.select(lambda x: x < new_age, ["age"]) self._evolve_particles(selection, selection.time_step + selection.age) return end_time = self.particles.time_step + self.particles.age - + self._evolve_particles(self.particles, end_time) - - + def commit_parameters(self): self.parameters.send_cached_parameters_to_code() self.overridden().commit_parameters() - + def initialize_module_with_current_parameters(self): self.parameters.send_cached_parameters_to_code() - + def initialize_module_with_default_parameters(self): self.parameters.set_defaults() self.commit_parameters() @@ -455,5 +560,6 @@ def initialize_module_with_default_parameters(self): def update_time_steps(self): pass - -Sse = SSE +# Backwards compatibility +SSE = Sse +SSEInterface = SseInterface From d6f9811d271fdc003bb55a41ad3a82917049d120 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Aug 2024 09:26:33 +0200 Subject: [PATCH 2/6] Update Makefile to use standardised name --- src/amuse/community/sse/Makefile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/amuse/community/sse/Makefile b/src/amuse/community/sse/Makefile index dc8d6ffa16..f901227d73 100644 --- a/src/amuse/community/sse/Makefile +++ b/src/amuse/community/sse/Makefile @@ -40,7 +40,7 @@ sse_worker: worker_code.f90 interface.o $(OBJ) $(MPIF90) $(FCFLAGS) $(FS_FLAGS) $(LDFLAGS) $^ -o $@ $(FS_LIBS) $(LIBS) worker_code.f90: interface.py - $(CODE_GENERATOR) --type=f90 interface.py SSEInterface -o $@ + $(CODE_GENERATOR) --type=f90 interface.py SseInterface -o $@ .f.o: $< $(FORTRAN) -c $(F77FLAGS) $(FCFLAGS) -o $@ $< From f0adf67a54f22b643ffcc89f095b1a22210b8ccc Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Aug 2024 10:15:25 +0200 Subject: [PATCH 3/6] update BSE interface for style --- src/amuse/community/bse/interface.py | 796 ++++++++++++++++----------- 1 file changed, 471 insertions(+), 325 deletions(-) diff --git a/src/amuse/community/bse/interface.py b/src/amuse/community/bse/interface.py index 4edad8aa17..884e14d87e 100644 --- a/src/amuse/community/bse/interface.py +++ b/src/amuse/community/bse/interface.py @@ -1,153 +1,286 @@ -from amuse.community import * -from amuse.units import units -from amuse.units import constants -from amuse.units.quantities import Quantity +""" +AMUSE interface for the binary population synthesis code BSE +""" + +import numpy + +from amuse.units import units, constants + from amuse.community.interface import common -from amuse.datamodel import Particles -from amuse.datamodel import ParticlesSubset +from amuse.datamodel import Particles, ParticlesSubset + +from amuse.support.interface import InCodeComponentImplementation +from amuse.community import ( + CodeInterface, + LiteratureReferencesMixIn, + legacy_function, + LegacyFunctionSpecification, + NO_UNIT, +) -import numpy -class BSEInterface(CodeInterface, common.CommonCodeInterface , LiteratureReferencesMixIn): +class BseInterface( + CodeInterface, common.CommonCodeInterface, LiteratureReferencesMixIn +): """ - Binary evolution is performed by the **rapid** binary-star evolution (BSE) - algorithm. Circularization of eccentric orbits and synchronization of stellar - rotation with the orbital motion owing to tidal interaction is modelled in detail. - Angular momentum loss mechanisms, such as gravitational radiation and magnetic - braking, are also modelled. Wind accretion, where the secondary may accrete some - of the material lost from the primary in a wind, is allowed with the necessary - adjustments made to the orbital parameters in the event of any mass variations. - Mass transfer also occurs if either star fills its Roche lobe and may proceed on a - nuclear, thermal or dynamical time-scale. In the latter regime, the radius of the - primary increases in response to mass-loss at a faster rate than the Roche-lobe of - the star. Stars with deep surface convection zones and degenerate stars are - unstable to such dynamical time-scale mass loss unless the mass ratio of the system - is less than some critical value. The outcome is a common-envelope event if the - primary is a giant star. This results in merging or formation of a close binary, or - a direct merging if the primary is a white dwarf or low-mass main-sequence star. On - the other hand, mass transfer on a nuclear or thermal time-scale is assumed to be a - steady process. Prescriptions to determine the type and rate of mass transfer, the - response of the secondary to accretion and the outcome of any merger events are in + Binary evolution is performed by the **rapid** binary-star evolution (BSE) + algorithm. Circularization of eccentric orbits and synchronization of stellar + rotation with the orbital motion owing to tidal interaction is modelled in detail. + Angular momentum loss mechanisms, such as gravitational radiation and magnetic + braking, are also modelled. Wind accretion, where the secondary may accrete some + of the material lost from the primary in a wind, is allowed with the necessary + adjustments made to the orbital parameters in the event of any mass variations. + Mass transfer also occurs if either star fills its Roche lobe and may proceed on a + nuclear, thermal or dynamical time-scale. In the latter regime, the radius of the + primary increases in response to mass-loss at a faster rate than the Roche-lobe of + the star. Stars with deep surface convection zones and degenerate stars are + unstable to such dynamical time-scale mass loss unless the mass ratio of the system + is less than some critical value. The outcome is a common-envelope event if the + primary is a giant star. This results in merging or formation of a close binary, or + a direct merging if the primary is a white dwarf or low-mass main-sequence star. On + the other hand, mass transfer on a nuclear or thermal time-scale is assumed to be a + steady process. Prescriptions to determine the type and rate of mass transfer, the + response of the secondary to accretion and the outcome of any merger events are in place in BSE and the details can be found in the BSE paper: - + .. [#] ADS:2002MNRAS.329..897H (Hurley J.R., Tout C.A., & Pols O.R., 2002, MNRAS, 329, 897) .. [#] ADS:2000MNRAS.315..543H (Hurley J.R., Pols O.R., Tout C.A., 2000, MNRAS, 315, 543) """ + def __init__(self, **options): CodeInterface.__init__(self, name_of_the_worker="bse_worker", **options) LiteratureReferencesMixIn.__init__(self) - @legacy_function + @legacy_function def initialize(): - function = LegacyFunctionSpecification() - function.addParameter('z_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('neta_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('bwind_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('hewind_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('alpha1_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('CElambda_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('ceflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('tflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('ifflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('wdflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('bhflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('nsflag_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('mxns_in', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('idum_in', dtype='i', direction=function.IN, unit = NO_UNIT) - function.addParameter('pts1_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('pts2_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('pts3_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('sigma_in', dtype='d', direction=function.IN, unit = units.km / units.s) - function.addParameter('beta_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('xi_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('acc2_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('epsnov_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('eddfac_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('gamma_in', dtype='d', direction=function.IN, unit = NO_UNIT) - function.addParameter('status', dtype='i', direction=function.OUT, unit = NO_UNIT) + function = LegacyFunctionSpecification() + function.addParameter("z_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("neta_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter( + "bwind_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "hewind_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "alpha1_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "CElambda_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "ceflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "tflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "ifflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "wdflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "bhflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "nsflag_in", dtype="i", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "mxns_in", dtype="d", direction=function.IN, unit=units.MSun + ) + function.addParameter("idum_in", dtype="i", direction=function.IN, unit=NO_UNIT) + function.addParameter("pts1_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("pts2_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("pts3_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter( + "sigma_in", dtype="d", direction=function.IN, unit=units.km / units.s + ) + function.addParameter("beta_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("xi_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter("acc2_in", dtype="d", direction=function.IN, unit=NO_UNIT) + function.addParameter( + "epsnov_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "eddfac_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter( + "gamma_in", dtype="d", direction=function.IN, unit=NO_UNIT + ) + function.addParameter("status", dtype="i", direction=function.OUT, unit=NO_UNIT) return function - - @legacy_function + + @legacy_function def evolve_binary(): function = LegacyFunctionSpecification() - function.can_handle_array = True - function.addParameter('type1', dtype='i', direction=function.INOUT, unit = units.stellar_type) - function.addParameter('type2', dtype='i', direction=function.INOUT, unit = units.stellar_type) - function.addParameter('initial_mass1', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('initial_mass2', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('mass1', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('mass2', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('radius1', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('radius2', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('luminosity1', dtype='d', direction=function.INOUT, unit = units.LSun) - function.addParameter('luminosity2', dtype='d', direction=function.INOUT, unit = units.LSun) - function.addParameter('core_mass1', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('core_mass2', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('core_radius1', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('core_radius2', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('convective_envelope_mass1', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('convective_envelope_mass2', dtype='d', direction=function.INOUT, unit = units.MSun) - function.addParameter('convective_envelope_radius1', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('convective_envelope_radius2', dtype='d', direction=function.INOUT, unit = units.RSun) - function.addParameter('spin1', dtype='d', direction=function.INOUT, unit = NO_UNIT) - function.addParameter('spin2', dtype='d', direction=function.INOUT, unit = NO_UNIT) - function.addParameter('epoch1', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('epoch2', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('MS_lifetime1', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('MS_lifetime2', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('age', dtype='d', direction=function.INOUT, unit = units.Myr) - function.addParameter('orbital_period', dtype='d', direction=function.INOUT, unit = units.day) - function.addParameter('eccentricity', dtype='d', direction=function.INOUT, unit = NO_UNIT) - function.addParameter('end_time', dtype='d', direction=function.INOUT, unit = units.Myr) + function.can_handle_array = True + function.addParameter( + "type1", dtype="i", direction=function.INOUT, unit=units.stellar_type + ) + function.addParameter( + "type2", dtype="i", direction=function.INOUT, unit=units.stellar_type + ) + function.addParameter( + "initial_mass1", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "initial_mass2", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "mass1", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "mass2", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "radius1", dtype="d", direction=function.INOUT, unit=units.RSun + ) + function.addParameter( + "radius2", dtype="d", direction=function.INOUT, unit=units.RSun + ) + function.addParameter( + "luminosity1", dtype="d", direction=function.INOUT, unit=units.LSun + ) + function.addParameter( + "luminosity2", dtype="d", direction=function.INOUT, unit=units.LSun + ) + function.addParameter( + "core_mass1", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "core_mass2", dtype="d", direction=function.INOUT, unit=units.MSun + ) + function.addParameter( + "core_radius1", dtype="d", direction=function.INOUT, unit=units.RSun + ) + function.addParameter( + "core_radius2", dtype="d", direction=function.INOUT, unit=units.RSun + ) + function.addParameter( + "convective_envelope_mass1", + dtype="d", + direction=function.INOUT, + unit=units.MSun, + ) + function.addParameter( + "convective_envelope_mass2", + dtype="d", + direction=function.INOUT, + unit=units.MSun, + ) + function.addParameter( + "convective_envelope_radius1", + dtype="d", + direction=function.INOUT, + unit=units.RSun, + ) + function.addParameter( + "convective_envelope_radius2", + dtype="d", + direction=function.INOUT, + unit=units.RSun, + ) + function.addParameter( + "spin1", dtype="d", direction=function.INOUT, unit=NO_UNIT + ) + function.addParameter( + "spin2", dtype="d", direction=function.INOUT, unit=NO_UNIT + ) + function.addParameter( + "epoch1", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter( + "epoch2", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter( + "MS_lifetime1", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter( + "MS_lifetime2", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter( + "age", dtype="d", direction=function.INOUT, unit=units.Myr + ) + function.addParameter( + "orbital_period", dtype="d", direction=function.INOUT, unit=units.day + ) + function.addParameter( + "eccentricity", dtype="d", direction=function.INOUT, unit=NO_UNIT + ) + function.addParameter( + "end_time", dtype="d", direction=function.INOUT, unit=units.Myr + ) return function - - @legacy_function + + @legacy_function def get_time_step(): - function = LegacyFunctionSpecification() + function = LegacyFunctionSpecification() function.can_handle_array = True - function.addParameter('type1', dtype='i', direction=function.IN, unit = units.stellar_type) - function.addParameter('type2', dtype='i', direction=function.IN, unit = units.stellar_type) - function.addParameter('initial_mass1', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('initial_mass2', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('mass1', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('mass2', dtype='d', direction=function.IN, unit = units.MSun) - function.addParameter('MS_lifetime1', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('MS_lifetime2', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('epoch1', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('epoch2', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('age', dtype='d', direction=function.IN, unit = units.Myr) - function.addParameter('time_step', dtype='d', direction=function.OUT, unit = units.Myr) - + function.addParameter( + "type1", dtype="i", direction=function.IN, unit=units.stellar_type + ) + function.addParameter( + "type2", dtype="i", direction=function.IN, unit=units.stellar_type + ) + function.addParameter( + "initial_mass1", dtype="d", direction=function.IN, unit=units.MSun + ) + function.addParameter( + "initial_mass2", dtype="d", direction=function.IN, unit=units.MSun + ) + function.addParameter( + "mass1", dtype="d", direction=function.IN, unit=units.MSun + ) + function.addParameter( + "mass2", dtype="d", direction=function.IN, unit=units.MSun + ) + function.addParameter( + "MS_lifetime1", dtype="d", direction=function.IN, unit=units.Myr + ) + function.addParameter( + "MS_lifetime2", dtype="d", direction=function.IN, unit=units.Myr + ) + function.addParameter( + "epoch1", dtype="d", direction=function.IN, unit=units.Myr + ) + function.addParameter( + "epoch2", dtype="d", direction=function.IN, unit=units.Myr + ) + function.addParameter("age", dtype="d", direction=function.IN, unit=units.Myr) + function.addParameter( + "time_step", dtype="d", direction=function.OUT, unit=units.Myr + ) + return function - + def get_time_step_for_binary(self, binary): - + current_values = {} - current_values['type1'] = binary.type1.value_in(units.stellar_type) - current_values['type2'] = binary.type2.value_in(units.stellar_type) - current_values['initial_mass1'] = binary.initial_mass1.value_in(units.MSun) - current_values['initial_mass2'] = binary.initial_mass2.value_in(units.MSun) - current_values['mass1'] = binary.mass1.value_in(units.MSun) - current_values['mass2'] = binary.mass2.value_in(units.MSun) - current_values['MS_lifetime1'] = binary.MS_lifetime1.value_in(units.Myr) - current_values['MS_lifetime2'] = binary.MS_lifetime2.value_in(units.Myr) - current_values['epoch1'] = binary.epoch1.value_in(units.Myr) - current_values['epoch2'] = binary.epoch2.value_in(units.Myr) - current_values['age'] = binary.age.value_in(units.Myr) - + current_values["type1"] = binary.type1.value_in(units.stellar_type) + current_values["type2"] = binary.type2.value_in(units.stellar_type) + current_values["initial_mass1"] = binary.initial_mass1.value_in(units.MSun) + current_values["initial_mass2"] = binary.initial_mass2.value_in(units.MSun) + current_values["mass1"] = binary.mass1.value_in(units.MSun) + current_values["mass2"] = binary.mass2.value_in(units.MSun) + current_values["MS_lifetime1"] = binary.MS_lifetime1.value_in(units.Myr) + current_values["MS_lifetime2"] = binary.MS_lifetime2.value_in(units.Myr) + current_values["epoch1"] = binary.epoch1.value_in(units.Myr) + current_values["epoch2"] = binary.epoch2.value_in(units.Myr) + current_values["age"] = binary.age.value_in(units.Myr) + result = self.get_time_step(**current_values) - + return result | units.Myr - - + def evolve_particle(self, particle, time_end): t = particle.current_time if particle.stellar_type == 15: return while t < time_end: t0 = t - t = t0 + self.get_time_step_for_binary(particle) + t = t0 + self.get_time_step_for_binary(particle) if t > time_end: t = time_end self.evolve_star(particle, t) @@ -155,52 +288,56 @@ def evolve_particle(self, particle, time_end): dt = t1 - t0 t0 = t1 if dt.value_in(units.Myr) == 0.0: - #print t, t0, t1, dt, "BREAK BREAK BREAK!" + # print t, t0, t1, dt, "BREAK BREAK BREAK!" return if particle.stellar_type == 15: return - + def initialize_code(self): return 0 - + def commit_parameters(self): return 0 - + def recommit_parameters(self): return 0 - + def cleanup_code(self): return 0 - + def commit_particles(self): return 0 - - - - -class BSEStars(Particles): - - def __init__(self, code_interface, storage = None): - Particles.__init__(self, storage = storage) - self._private.code_interface = code_interface - self.add_calculated_attribute("temperature", self.calculate_effective_temperature, ["luminosity", "radius"]) - + + +class BseStars(Particles): + + def __init__(self, code_interface, storage=None): + Particles.__init__(self, storage=storage) + self._private.code_interface = code_interface + self.add_calculated_attribute( + "temperature", + self.calculate_effective_temperature, + ["luminosity", "radius"], + ) + def calculate_effective_temperature(self, luminosity, radius): - return ((luminosity/(constants.four_pi_stefan_boltzmann*radius**2))**.25).in_(units.K) - - def add_particles_to_store(self, keys, attributes = [], values = []): + return ( + (luminosity / (constants.four_pi_stefan_boltzmann * radius**2)) ** 0.25 + ).in_(units.K) + + def add_particles_to_store(self, keys, attributes=[], values=[]): if len(keys) == 0: return - + all_attributes = [] all_attributes.extend(attributes) all_values = [] all_values.extend(values) - + mapping_from_attribute_to_default_value = { - "stellar_type" : 1 | units.stellar_type, + "stellar_type": 1 | units.stellar_type, "radius": 0 | units.RSun, - "luminosity": 0 | units.LSun, + "luminosity": 0 | units.LSun, "core_mass": 0 | units.MSun, "core_radius": 0 | units.RSun, "convective_envelope_mass": 0 | units.MSun, @@ -209,370 +346,375 @@ def add_particles_to_store(self, keys, attributes = [], values = []): "spin": 0 | units.none, "main_sequence_lifetime": 0 | units.Myr, "age": 0 | units.Myr, - "stellar_type": 0 | units.stellar_type #units.stellar_type("Main Sequence star"), + "stellar_type": 0 + | units.stellar_type, # units.stellar_type("Main Sequence star"), } - + given_attributes = set(attributes) - + if not "initial_mass" in given_attributes: index_of_mass_attibute = attributes.index("mass") all_attributes.append("initial_mass") all_values.append(values[index_of_mass_attibute] * 1.0) - + for attribute, default_value in mapping_from_attribute_to_default_value.items(): if not attribute in given_attributes: all_attributes.append(attribute) all_values.append(default_value.as_vector_with_length(len(keys))) - - super(BSEStars, self).add_particles_to_store(keys, all_attributes, all_values) - - + + super(BseStars, self).add_particles_to_store(keys, all_attributes, all_values) + def get_defined_attribute_names(self): return ["mass", "radius"] -class BSEBinaries(Particles): - - def __init__(self, code_interface, storage = None): - Particles.__init__(self, storage = storage) - self._private.code_interface = code_interface - - def add_particles_to_store(self, keys, attributes = [], values = []): + +class BseBinaries(Particles): + + def __init__(self, code_interface, storage=None): + Particles.__init__(self, storage=storage) + self._private.code_interface = code_interface + + def add_particles_to_store(self, keys, attributes=[], values=[]): if len(keys) == 0: return - + given_attributes = set(attributes) - + if not "child1" in given_attributes: raise Exception("a binary must always have a child1 attribute") - + if not "child2" in given_attributes: raise Exception("a binary must always have a child2 attribute") - + all_attributes = [] all_values = [] for attribute, value in zip(attributes, values): all_attributes.append(attribute) - if attribute == 'child1' or attribute == 'child2': - value = value.copy_with_link_transfer(None, self._private.code_interface.particles) + if attribute == "child1" or attribute == "child2": + value = value.copy_with_link_transfer( + None, self._private.code_interface.particles + ) all_values.append(value) else: all_values.append(value) - + mapping_from_attribute_to_default_value = { "eccentricity": 0.0 | units.none, - "age": 0 | units.Myr + "age": 0 | units.Myr, } - - - + for attribute, default_value in mapping_from_attribute_to_default_value.items(): if not attribute in given_attributes: all_attributes.append(attribute) all_values.append(default_value.as_vector_with_length(len(keys))) - - super(BSEBinaries, self).add_particles_to_store(keys, all_attributes, all_values) - + + super(BseBinaries, self).add_particles_to_store( + keys, all_attributes, all_values + ) + added_particles = ParticlesSubset(self, keys) self._private.code_interface._evolve_binaries(added_particles, 1e-08 | units.yr) - + def get_defined_attribute_names(self): return ["eccentricity", "orbital_period", "age", "child1", "child2"] - -class BSE(common.CommonCode): - + + +class Bse(common.CommonCode): + def __init__(self, **options): - InCodeComponentImplementation.__init__(self, BSEInterface(**options), **options) + InCodeComponentImplementation.__init__(self, BseInterface(**options), **options) self.model_time = 0.0 | units.yr - - + def define_parameters(self, handler): - + handler.add_caching_parameter( - "initialize", - "z_in", - "metallicity", - "Metallicity of all stars", - 0.02 + "initialize", "z_in", "metallicity", "Metallicity of all stars", 0.02 ) - + handler.add_caching_parameter( "initialize", "neta_in", "reimers_mass_loss_coefficient", "Reimers mass-loss coefficient (neta*4x10^-13; 0.5 normally)", - 0.5 + 0.5, ) - + handler.add_caching_parameter( "initialize", "bwind_in", "binary_enhanced_mass_loss_parameter", "The binary enhanced mass loss parameter (inactive for single).", - 0.0 + 0.0, ) - + handler.add_caching_parameter( "initialize", "hewind_in", "helium_star_mass_loss_factor", "Helium star mass loss factor", - 1.0 + 1.0, ) - + handler.add_caching_parameter( "initialize", "alpha1_in", "common_envelope_efficiency", "The common-envelope efficiency parameter", - 1.0 + 1.0, ) - + handler.add_caching_parameter( "initialize", "CElambda_in", "common_envelope_binding_energy_factor", "The binding energy factor for common envelope evolution", - 0.5 + 0.5, ) - + handler.add_caching_parameter( "initialize", "ceflag_in", "common_envelope_model_flag", "ceflag > 0 activates spin-energy correction in common-envelope. ceflag = 3 activates de Kool common-envelope model (0).", - 0 + 0, ) - + handler.add_caching_parameter( "initialize", "tflag_in", "tidal_circularisation_flag", "tflag > 0 activates tidal circularisation (1).", - 1 + 1, ) - + handler.add_caching_parameter( "initialize", "ifflag_in", "white_dwarf_IFMR_flag", "ifflag > 0 uses white dwarf IFMR (initial-final mass relation) of HPE, 1995, MNRAS, 272, 800 (0).", - 0 + 0, ) - + handler.add_caching_parameter( "initialize", "wdflag_in", "white_dwarf_cooling_flag", "wdflag > 0 uses modified-Mestel cooling for WDs (0).", - 1 + 1, ) - + handler.add_caching_parameter( "initialize", "bhflag_in", "black_hole_kick_flag", "bhflag > 0 allows velocity kick at BH formation (0).", - 0 + 0, ) - + handler.add_caching_parameter( "initialize", "nsflag_in", "neutron_star_mass_flag", "nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1).", - 1 + 1, ) - + handler.add_caching_parameter( "initialize", "mxns_in", "maximum_neutron_star_mass", "The maximum neutron star mass (1.8, nsflag=0; 3.0, nsflag=1).", - 3.0 | units.MSun + 3.0 | units.MSun, ) - + handler.add_caching_parameter( "initialize", "idum_in", "SN_kick_random_seed", "The random number seed used in the kick routine.", - 29769 + 29769, ) - + handler.add_caching_parameter( "initialize", "pts1_in", "fractional_time_step_1", "The timesteps chosen in each evolution phase as decimal fractions of the time taken in that phase: MS (0.05)", - 0.05 + 0.05, ) - + handler.add_caching_parameter( "initialize", "pts2_in", "fractional_time_step_2", "The timesteps chosen in each evolution phase as decimal fractions of the time taken in that phase: GB, CHeB, AGB, HeGB (0.01)", - 0.01 + 0.01, ) - + handler.add_caching_parameter( "initialize", "pts3_in", "fractional_time_step_3", "The timesteps chosen in each evolution phase as decimal fractions of the time taken in that phase: HG, HeMS (0.02)", - 0.02 + 0.02, ) - + handler.add_caching_parameter( "initialize", "sigma_in", "SN_kick_speed_dispersion", "The dispersion in the Maxwellian for the SN kick speed (190 km/s).", - 190.0 | units.km / units.s + 190.0 | units.km / units.s, ) - + handler.add_caching_parameter( "initialize", "beta_in", "wind_velocity_factor", "The wind velocity factor: proportional to vwind**2 (1/8).", - 0.125 + 0.125, ) - + handler.add_caching_parameter( "initialize", "xi_in", "wind_accretion_efficiency", "The wind accretion efficiency factor (1.0).", - 1.0 + 1.0, ) - + handler.add_caching_parameter( "initialize", "acc2_in", "wind_accretion_factor", "The Bondi-Hoyle wind accretion factor (3/2).", - 1.5 + 1.5, ) - + handler.add_caching_parameter( "initialize", "epsnov_in", "nova_retained_accreted_matter_fraction", "The fraction of accreted matter retained in nova eruption (0.001).", - 0.001 + 0.001, ) - + handler.add_caching_parameter( "initialize", "eddfac_in", "Eddington_mass_transfer_limit_factor", "The Eddington limit factor for mass transfer (1.0).", - 1.0 + 1.0, ) - + handler.add_caching_parameter( "initialize", "gamma_in", "Roche_angular_momentum_factor", "The angular momentum factor for mass lost during Roche (-1.0). ", - -1.0 + -1.0, ) - - + def define_state(self, handler): common.CommonCode.define_state(self, handler) - handler.add_transition('INITIALIZED','RUN','commit_parameters') - handler.add_method('RUN', 'evolve_binary') - - handler.add_method('RUN','before_get_parameter') - handler.add_method('RUN','before_set_parameter') - - - - + handler.add_transition("INITIALIZED", "RUN", "commit_parameters") + handler.add_method("RUN", "evolve_binary") + + handler.add_method("RUN", "before_get_parameter") + handler.add_method("RUN", "before_set_parameter") + def define_particle_sets(self, handler): - handler.define_inmemory_set('particles', BSEStars) - handler.define_inmemory_set('binaries', BSEBinaries) - + handler.define_inmemory_set("particles", BseStars) + handler.define_inmemory_set("binaries", BseBinaries) + handler.add_attribute( - 'binaries', - 'time_step', - '_get_time_step', - ('child1', 'child2', 'age') - #('child1', 'type2', - # 'initial_mass1', 'initial_mass2', + "binaries", + "time_step", + "_get_time_step", + ("child1", "child2", "age"), + # ('child1', 'type2', + # 'initial_mass1', 'initial_mass2', # 'mass1', 'mass2', # 'MS_lifetime1', 'MS_lifetime2', # 'epoch1', 'epoch2', #'age') ) - + def _get_time_step(self, child1, child2, age): child1 = child1.as_set() child2 = child2.as_set() return self.get_time_step( - child1.stellar_type, child2.stellar_type, - child1.initial_mass, child2.initial_mass, - child1.mass, child2.mass, - child1.age, child2.age, - child1.epoch, child2.epoch, - age - ) - + child1.stellar_type, + child2.stellar_type, + child1.initial_mass, + child2.initial_mass, + child1.mass, + child2.mass, + child1.age, + child2.age, + child1.epoch, + child2.epoch, + age, + ) + def orbital_period_to_semi_major_axis(self, orbital_period, mass1, mass2): mu = (mass1 + mass2) * constants.G - return (((orbital_period / (2.0 * numpy.pi))**2)*mu)**(1.0/3.0) - + return (((orbital_period / (2.0 * numpy.pi)) ** 2) * mu) ** (1.0 / 3.0) + def semi_major_axis_to_orbital_period(self, semi_major_axis, mass1, mass2): mu = (mass1 + mass2) * constants.G - return 2.0 * numpy.pi * ((semi_major_axis**3/mu)**0.5) - + return 2.0 * numpy.pi * ((semi_major_axis**3 / mu) ** 0.5) + def _evolve_binaries(self, particles, end_time): - binary_attributes = ( - "age", - "semi_major_axis", - "eccentricity" - ) - + binary_attributes = ("age", "semi_major_axis", "eccentricity") + single_attributes = ( "stellar_type", - "initial_mass", - "mass", - "radius", + "initial_mass", + "mass", + "radius", "luminosity", - "core_mass", - "core_radius", - "convective_envelope_mass", - "convective_envelope_radius", - "spin", - "epoch", - "age", - ) - + "core_mass", + "core_radius", + "convective_envelope_mass", + "convective_envelope_radius", + "spin", + "epoch", + "age", + ) + children1 = particles.child1.as_set() children2 = particles.child2.as_set() - children1_arguments = children1.get_values_in_store(children1.get_all_indices_in_store(), single_attributes) - children2_arguments = children2.get_values_in_store(children2.get_all_indices_in_store(), single_attributes) - - binaries_arguments = particles.get_values_in_store(particles.get_all_indices_in_store(), binary_attributes) - - binaries_arguments[1] = self.semi_major_axis_to_orbital_period(binaries_arguments[1] , children1_arguments[2], children2_arguments[2]) - + children1_arguments = children1.get_values_in_store( + children1.get_all_indices_in_store(), single_attributes + ) + children2_arguments = children2.get_values_in_store( + children2.get_all_indices_in_store(), single_attributes + ) + + binaries_arguments = particles.get_values_in_store( + particles.get_all_indices_in_store(), binary_attributes + ) + + binaries_arguments[1] = self.semi_major_axis_to_orbital_period( + binaries_arguments[1], children1_arguments[2], children2_arguments[2] + ) + arguments = [] for argument1, argument2 in zip(children1_arguments, children2_arguments): arguments.append(argument1) arguments.append(argument2) - + arguments.extend(binaries_arguments) arguments.append(end_time.as_vector_with_length(len(particles))) - + result = self.evolve_binary(*arguments) - - result[-3] = self.orbital_period_to_semi_major_axis(result[-3] , result[4], result[5]) - - + + result[-3] = self.orbital_period_to_semi_major_axis( + result[-3], result[4], result[5] + ) + children1_results = [] children2_results = [] index = 0 @@ -581,48 +723,57 @@ def _evolve_binaries(self, particles, end_time): index += 1 children2_results.append(result[index]) index += 1 - - children1.set_values_in_store(children1.get_all_indices_in_store(), single_attributes, children1_results) - children2.set_values_in_store(children2.get_all_indices_in_store(), single_attributes, children2_results) - particles.set_values_in_store(particles.get_all_indices_in_store(), binary_attributes, result[index:]) - - - def evolve_model(self, end_time = None, keep_synchronous = True): + + children1.set_values_in_store( + children1.get_all_indices_in_store(), single_attributes, children1_results + ) + children2.set_values_in_store( + children2.get_all_indices_in_store(), single_attributes, children2_results + ) + particles.set_values_in_store( + particles.get_all_indices_in_store(), binary_attributes, result[index:] + ) + + def evolve_model(self, end_time=None, keep_synchronous=True): if not keep_synchronous: - self._evolve_binaries(self.binaries, self.binaries.time_step + self.binaries.age) + self._evolve_binaries( + self.binaries, self.binaries.time_step + self.binaries.age + ) return - + if end_time is None: end_time = self.model_time + min(self.binaries.time_step) - self._evolve_binaries(self.binaries, end_time - self.model_time + self.binaries.age) + self._evolve_binaries( + self.binaries, end_time - self.model_time + self.binaries.age + ) self.model_time = end_time - + def commit_particles(self): pass def update_time_steps(self): pass - + def commit_parameters(self): self.parameters.send_cached_parameters_to_code() self.overridden().commit_parameters() - + def initialize_module_with_current_parameters(self): self.commit_parameters() - + def initialize_module_with_default_parameters(self): """ * neta is the Reimers mass-loss coefficent (neta*4x10^-13; 0.5 normally). * bwind is the binary enhanced mass loss parameter (inactive for single). - * hewind is a helium star mass loss factor (1.0 normally). - * sigma is the dispersion in the Maxwellian for the SN kick speed (190 km/s). + * hewind is a helium star mass loss factor (1.0 normally). + * sigma is the dispersion in the Maxwellian for the SN kick speed (190 km/s). * - * ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0). - * wdflag > 0 uses modified-Mestel cooling for WDs (0). - * bhflag > 0 allows velocity kick at BH formation (0). - * nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1). - * mxns is the maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1). - * idum is the random number seed used in the kick routine. + * ifflag > 0 uses WD IFMR of HPE, 1995, MNRAS, 272, 800 (0). + * wdflag > 0 uses modified-Mestel cooling for WDs (0). + * bhflag > 0 allows velocity kick at BH formation (0). + * nsflag > 0 takes NS/BH mass from Belczynski et al. 2002, ApJ, 572, 407 (1). + * mxns is the maximum NS mass (1.8, nsflag=0; 3.0, nsflag=1). + * idum is the random number seed used in the kick routine. * * Next come the parameters that determine the timesteps chosen in each * evolution phase: @@ -634,11 +785,6 @@ def initialize_module_with_default_parameters(self): self.parameters.set_defaults() self.commit_parameters() - - - - - - - -Bse = BSE +# backwards compatibility +BSEInterface = BseInterface +BSE = Bse From 7b4f8190989e87b16d3a25255344d9a35a34e49e Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Aug 2024 10:20:14 +0200 Subject: [PATCH 4/6] Cleanup fortran interface, remove old code --- src/amuse/community/bse/Makefile | 2 +- src/amuse/community/bse/interface.f | 15 ++------------- 2 files changed, 3 insertions(+), 14 deletions(-) diff --git a/src/amuse/community/bse/Makefile b/src/amuse/community/bse/Makefile index 72e3aa454d..d85e7878c4 100644 --- a/src/amuse/community/bse/Makefile +++ b/src/amuse/community/bse/Makefile @@ -40,7 +40,7 @@ bse_worker: worker_code.f90 interface.o $(BSEOBJ) $(MPIFC) $(FFLAGS) $(FS_FLAGS) $(LDFLAGS) $^ -o $@ $(FS_LIBS) $(LIBS) worker_code.f90: interface.py - $(CODE_GENERATOR) --type=f90 interface.py BSEInterface -o $@ + $(CODE_GENERATOR) --type=f90 interface.py BseInterface -o $@ .f.o: $< $(FORTRAN) -c $(F77FLAGS) $(FFLAGS) -o $@ $< diff --git a/src/amuse/community/bse/interface.f b/src/amuse/community/bse/interface.f index 31e10a132b..352d297f07 100644 --- a/src/amuse/community/bse/interface.f +++ b/src/amuse/community/bse/interface.f @@ -20,7 +20,7 @@ subroutine initialize(z_in, real*8 mxns_in, pts1_in, pts2_in, pts3_in real*8 sigma_in, beta_in, xi_in, acc2_in real*8 epsnov_in, eddfac_in, gamma_in - CHARACTER*8 label(14) + character*8 label(14) integer status include 'src/const_bse.h' common /SSE_init/ z, zpars @@ -61,7 +61,7 @@ subroutine initialize(z_in, * * Set the collision matrix. * - CALL instar + call instar * label(1) = 'INITIAL ' label(2) = 'KW CHNGE' @@ -145,15 +145,6 @@ subroutine evolve_binary(type1,type2,initial_mass1,initial_mass2, ! dtp = 0.0d0 ! dtp = age+1 -c call evolv2(type1,type2,initial_mass1,initial_mass2, -c & mass1, mass2, radius1, radius2, luminosity1, luminosity2, -c & core_mass1, core_mass2, core_radius1, core_radius2, -c & envelope_mass1, envelope_mass2, envelope_radius1, -c & envelope_radius2, spin1, spin2, epoch1, epoch2, -c & MS_lifetime1, MS_lifetime2, age, end_time, -c & dtp,z,zpars, -c & orbital_period, eccentricity) - call evolv2(kstar,mass0,mass,rad,lum,massc,radc, & menv,renv,ospin,epoch,tms, & age, end_time,dtp,z,zpars, @@ -193,8 +184,6 @@ subroutine get_time_step(type1, type2, initial_mass1, & initial_mass2, mass1, mass2, MS_lifetime1, & MS_lifetime2, epoch1, epoch2, age, time_step) -cf2py intent(out) dt -cf2py intent(in) kw, mass, age, mt, tm, epoch implicit none integer type1, type2 real*8 initial_mass1, initial_mass2, mass1, mass2 From 1b6365f70f2fe864a46687510e92a8a6fe89d585 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Aug 2024 10:23:36 +0200 Subject: [PATCH 5/6] update sse/bse tests --- src/amuse/test/suite/codes_tests/test_bse.py | 36 +++++------ src/amuse/test/suite/codes_tests/test_sse.py | 64 ++++++++++---------- 2 files changed, 50 insertions(+), 50 deletions(-) diff --git a/src/amuse/test/suite/codes_tests/test_bse.py b/src/amuse/test/suite/codes_tests/test_bse.py index da464514e2..b5c6f28c5e 100644 --- a/src/amuse/test/suite/codes_tests/test_bse.py +++ b/src/amuse/test/suite/codes_tests/test_bse.py @@ -1,4 +1,4 @@ -from amuse.community.bse.interface import BSE, BSEInterface +from amuse.community.bse.interface import Bse, BseInterface from amuse.test.amusetest import TestWithMPI from amuse.units import units @@ -8,7 +8,7 @@ import numpy -class TestBSEInterface(TestWithMPI): +class TestBseInterface(TestWithMPI): class state(object): def __init__(self): @@ -42,7 +42,7 @@ def __init__(self): def test1(self): print("Test initialization...") - instance = BSEInterface() + instance = BseInterface() metallicity = 0.02 neta = 0.5 bwind = 0.0 @@ -78,7 +78,7 @@ def test1(self): def test2(self): print("Test basic operations (legacy functions evolve & get_time_step)...") - instance = BSEInterface() + instance = BseInterface() status = instance.initialize(0.02, 0.5, 0.0, 0.5, 1.0, 0.5, 0, 1, 0, 1, 0, 1, 3.0, 29769, 0.05, 0.01, 0.02, 190.0, 1.0/8.0, 1.0, 3.0/2.0, 0.001, 1.0, -1.0) @@ -145,7 +145,7 @@ def test2(self): def test3(self): print("Test whether the interface can handle arrays...") - instance = BSEInterface() + instance = BseInterface() status = instance.initialize(0.02, 0.5, 0.0, 0.5, 1.0, 0.5, 0, 1, 0, 1, 0, 1, 3.0, 29769, 0.05, 0.01, 0.02, 190.0, 1.0/8.0, 1.0, 3.0/2.0, 0.001, 1.0, -1.0) masses1 = [10.0, 5.0, 4.0] @@ -179,7 +179,7 @@ def test3(self): def test4(self): print("Test large number of particles...") number_of_particles = 2000 - instance = BSEInterface() + instance = BseInterface() status = instance.initialize(0.02, 0.5, 0.0, 0.5, 1.0, 0.5, 0, 1, 0, 1, 0, 1, 3.0, 29769, 0.05, 0.01, 0.02, 190.0, 1.0/8.0, 1.0, 3.0/2.0, 0.001, 1.0, -1.0) masses1 = [1.0 + ((x / 1.0*number_of_particles) * 10.0) for x in range(1, number_of_particles+1)] @@ -210,11 +210,11 @@ def test4(self): instance.stop() -class TestBSE(TestWithMPI): +class TestBse(TestWithMPI): def test1(self): print("Testing evolution of a close binary system...") - instance = BSE() + instance = Bse() instance.initialize_code() instance.parameters.metallicity = 0.001 instance.parameters.common_envelope_efficiency = 3.0 @@ -300,7 +300,7 @@ def test1(self): def test2(self): print("Testing evolution of a wide binary system.") - instance = BSE() + instance = Bse() instance.parameters.metallicity = 0.001 instance.parameters.common_envelope_efficiency = 3.0 instance.parameters.Eddington_mass_transfer_limit_factor = 10.0 @@ -385,7 +385,7 @@ def test2(self): def test3(self): print("Testing standard BSE example 2...") - instance = BSE() + instance = Bse() instance.parameters.common_envelope_efficiency = 3.0 instance.parameters.Eddington_mass_transfer_limit_factor = 10.0 instance.commit_parameters() @@ -477,7 +477,7 @@ def test3(self): def test4(self): print("Quick testing standard BSE example 2...") - instance = BSE() + instance = Bse() instance.parameters.common_envelope_efficiency = 3.0 instance.parameters.Eddington_mass_transfer_limit_factor = 10.0 instance.commit_parameters() @@ -519,7 +519,7 @@ def test4(self): def test5(self): print("Testing stellar collision...") - instance = BSE() + instance = Bse() instance.parameters.common_envelope_efficiency = 3.0 instance.parameters.Eddington_mass_transfer_limit_factor = 10.0 instance.commit_parameters() @@ -562,7 +562,7 @@ def test5(self): def test6(self): print("Testing additional parameters for initialization...") - instance = BSE() + instance = Bse() instance.initialize_code() self.assertEqual(instance.parameters.reimers_mass_loss_coefficient, 0.5) myvalue = 0.7 @@ -572,7 +572,7 @@ def test6(self): self.assertEqual(instance.parameters.reimers_mass_loss_coefficient, myvalue) instance.stop() - instance = BSE() + instance = Bse() self.assertEqual(instance.parameters.reimers_mass_loss_coefficient, 0.5) myvalue = 0.7 instance.parameters.reimers_mass_loss_coefficient = myvalue @@ -584,7 +584,7 @@ def test6(self): def test7(self): print("Test evolve_model optional arguments: end_time and keep_synchronous") - instance = BSE() + instance = Bse() instance.commit_parameters() stars = Particles(6) @@ -631,7 +631,7 @@ def test7(self): def test8(self): print("Testing adding and removing particles from stellar evolution code...") - instance = BSE() + instance = Bse() instance.initialize_code() stars = Particles(6) @@ -679,7 +679,7 @@ def test8(self): def test9(self): print("Testing BSE states") - instance = BSE() + instance = Bse() stars = Particles(2) stars.mass = [1.0, 0.2] | units.MSun @@ -705,7 +705,7 @@ def test9(self): print("initialize_code(), commit_parameters(), " "and cleanup_code() should be called automatically:", end=' ') - instance = BSE() + instance = Bse() self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') instance.parameters.reimers_mass_loss_coefficient = 0.5 self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') diff --git a/src/amuse/test/suite/codes_tests/test_sse.py b/src/amuse/test/suite/codes_tests/test_sse.py index 1d41fe9254..0fbde2b209 100644 --- a/src/amuse/test/suite/codes_tests/test_sse.py +++ b/src/amuse/test/suite/codes_tests/test_sse.py @@ -5,7 +5,7 @@ import numpy from subprocess import call -from amuse.community.sse.interface import SSEInterface, SSE +from amuse.community.sse.interface import SseInterface, Sse from amuse.test.amusetest import get_path_to_results, TestWithMPI from amuse import io @@ -56,7 +56,7 @@ def initialize_module_with_default_parameters(self, sse): pts1, pts2, pts3) def test1(self): - sse = SSEInterface() + sse = SseInterface() metallicity = 0.02 @@ -85,7 +85,7 @@ def test1(self): sse.stop() def test2(self): - sse = SSEInterface() + sse = SseInterface() metallicity = 0.02 @@ -159,7 +159,7 @@ def test2(self): sse.stop() def test3(self): - sse = SSEInterface() + sse = SseInterface() self.initialize_module_with_default_parameters(sse) types = [1, 1, 1] masses = [10, 5, 4] @@ -189,7 +189,7 @@ def test3(self): sse.stop() def test4(self): - sse = SSEInterface() + sse = SseInterface() self.initialize_module_with_default_parameters(sse) types = [1 for x in range(1, 4000)] masses = [1.0 + ((x / 4000.0) * 10.0) for x in range(1, 4000)] @@ -219,10 +219,10 @@ def test4(self): sse.stop() -class TestSSE(TestWithMPI): +class TestSse(TestWithMPI): def test1(self): - sse = SSE() + sse = Sse() sse.commit_parameters() stars = Particles(1) star = stars[0] @@ -290,7 +290,7 @@ def test1(self): sse.stop() def test2(self): - sse = SSE() + sse = Sse() sse.commit_parameters() stars = Particles(1) @@ -307,7 +307,7 @@ def test2(self): sse.stop() def test3(self): - sse = SSE() + sse = Sse() sse.commit_parameters() stars = Particles(1) @@ -332,7 +332,7 @@ def test3(self): sse.stop() def test5(self): - sse = SSE() + sse = Sse() sse.commit_parameters() stars = Particles(1) @@ -367,7 +367,7 @@ def test6(self): stars.mass = masses # Initialize stellar evolution code - instance = SSE() + instance = Sse() instance.commit_parameters() instance.particles.add_particles(stars) instance.commit_particles() @@ -398,7 +398,7 @@ def test7(self): stars = Particles(2) stars.mass = 1.0 | units.MSun for star in stars: - stellar_evolution = SSE() + stellar_evolution = Sse() stellar_evolution.commit_parameters() stellar_evolution.particles.add_particles(star.as_set()) stellar_evolution.commit_particles() @@ -412,7 +412,7 @@ def test7(self): print("Solved: SSE_muse_interface.f sets initial_mass to mass when necessary.") def test8(self): - instance = SSE() + instance = Sse() self.assertEqual(instance.parameters.reimers_mass_loss_coefficient, 0.5) myvalue = 0.7 instance.parameters.reimers_mass_loss_coefficient = myvalue @@ -421,7 +421,7 @@ def test8(self): self.assertEqual(instance.parameters.reimers_mass_loss_coefficient, myvalue) instance.stop() - instance = SSE() + instance = Sse() self.assertEqual(instance.parameters.reimers_mass_loss_coefficient, 0.5) myvalue = 0.7 instance.parameters.reimers_mass_loss_coefficient = myvalue @@ -432,7 +432,7 @@ def test8(self): def test9(self): print("Test: large number of particles") - stellar_evolution = SSE(max_message_length=500) + stellar_evolution = Sse(max_message_length=500) stellar_evolution.commit_parameters() number_of_particles = 10000 print("Has been tested with up to a million particles!") @@ -444,7 +444,7 @@ def test9(self): stellar_evolution.stop() def test10(self): - stellar_evolution = SSE() + stellar_evolution = Sse() stellar_evolution.commit_parameters() stars = Particles(10) stars.mass = 1.0 | units.MSun @@ -465,7 +465,7 @@ def test11(self): print("Test evolve_model optional arguments: end_time and keep_synchronous") stars = Particles(3) stars.mass = [1.0, 2.0, 3.0] | units.MSun - instance = SSE() + instance = Sse() instance.commit_parameters() instance.particles.add_particles(stars) @@ -498,7 +498,7 @@ def test12(self): particles = Particles(3) particles.mass = 1.0 | units.MSun - instance = SSE() + instance = Sse() instance.initialize_code() instance.commit_parameters() self.assertEqual(len(instance.particles), 0) # before creation @@ -531,7 +531,7 @@ def test13(self): stars.mass = 1.0 | units.MSun print("First do everything manually:", end=' ') - instance = SSE() + instance = Sse() self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') instance.initialize_code() self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') @@ -544,7 +544,7 @@ def test13(self): print("initialize_code(), commit_parameters(), " "and cleanup_code() should be called automatically:", end=' ') - instance = SSE() + instance = Sse() self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') instance.parameters.reimers_mass_loss_coefficient = 0.5 self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') @@ -558,7 +558,7 @@ def test14a(self): print("Testing basic operations: evolve_one_step and evolve_for (on particle)") stars = Particles(2) stars.mass = 1.0 | units.MSun - instance = SSE() + instance = Sse() se_stars = instance.particles.add_particles(stars) self.assertAlmostEqual(se_stars.age, [0.0, 0.0] | units.yr) @@ -580,7 +580,7 @@ def test14b(self): print("Testing basic operations: evolve_one_step and evolve_for (on subset)") stars = Particles(2) stars.mass = 1.0 | units.MSun - instance = SSE() + instance = Sse() se_stars = instance.particles.add_particles(stars) self.assertAlmostEqual(se_stars.age, [0.0, 0.0] | units.yr) @@ -619,7 +619,7 @@ def random_sample(self, N): stars = Particles(mass=masses) - instance = SSE() + instance = Sse() instance.particles.add_particles(stars) i = 0 @@ -649,7 +649,7 @@ def random_sample(self, N): stars = Particles(mass=masses) - instance = SSE() + instance = Sse() instance.particles.add_particles(stars) i = 0 @@ -662,7 +662,7 @@ def test17(self): print("evolve_one_step and evolve_for after particle removal and addition") particles = Particles(10) particles.mass = range(1, 11) | units.MSun - instance = SSE() + instance = Sse() instance.particles.add_particles(particles) self.assertAlmostEqual(instance.particles.age, 0.0 | units.yr) time_steps = numpy.linspace(0.1, 1.0, num=10) | units.Myr @@ -689,10 +689,10 @@ def test17(self): def test18(self): print("SSE validation") - sse_src_path = os.path.join(os.path.dirname(sys.modules[SSE.__module__].__file__), 'src') + sse_src_path = os.path.join(os.path.dirname(sys.modules[Sse.__module__].__file__), 'src') if not os.path.exists(os.path.join(sse_src_path, "evolve.in")): self.skip("Not in a source release") - instance = SSE() + instance = Sse() instance.particles.add_particle(Particle(mass=1.416 | units.MSun)) instance.particles[0].evolve_for(7000.0 | units.Myr) evolved_star = instance.particles.copy()[0] @@ -723,7 +723,7 @@ def test18(self): def test19(self): print("SSE core_mass and CO_core_mass (high mass star)") - instance = SSE() + instance = Sse() star = instance.particles.add_particle(Particle(mass=30 | units.MSun)) instance.evolve_model(5.8 | units.Myr) print(star.mass, star.core_mass, star.CO_core_mass, star.stellar_type) @@ -759,7 +759,7 @@ def test19(self): def test20(self): print("SSE core_mass and CO_core_mass (low mass stars)") - instance = SSE() + instance = Sse() stars = instance.particles.add_particles(Particles(mass=[0.6, 1.0] | units.MSun)) instance.evolve_model(100 | units.Gyr) self.assertEqual(str(stars[0].stellar_type), "Helium White Dwarf") @@ -773,7 +773,7 @@ def test20(self): instance.stop() def test21(self): - instance = SSE() + instance = Sse() stars = instance.particles.add_particles(Particles(mass=30 | units.MSun)) mass_loss_wind = stars[0].mass_loss_wind self.assertAlmostRelativeEquals(mass_loss_wind, 1.703e-07 | units.MSun / units.yr, 3) @@ -785,7 +785,7 @@ def test21(self): instance.stop() def test22(self): - instance = SSE() + instance = Sse() stars = instance.particles.add_particles(Particles(mass=[1.0, 10.0] | units.MSun)) gyration_radius = stars.gyration_radius self.assertTrue(numpy.all(0.0 < gyration_radius)) @@ -798,7 +798,7 @@ def test22(self): instance.stop() def test23(self): - instance = SSE() + instance = Sse() p = Particles(mass=[1.0, 10.0] | units.MSun, temperature=[10, 10] | units.K) stars = instance.particles.add_particles(p) channel = stars.new_channel_to(p) From 9f66ae027e7b66cf2ae8356053ac994b928dbab4 Mon Sep 17 00:00:00 2001 From: Steven Rieder Date: Thu, 22 Aug 2024 10:26:29 +0200 Subject: [PATCH 6/6] pass sse/bse tests through Black --- src/amuse/test/suite/codes_tests/test_bse.py | 574 ++++++++++++++----- src/amuse/test/suite/codes_tests/test_sse.py | 415 ++++++++++---- 2 files changed, 728 insertions(+), 261 deletions(-) diff --git a/src/amuse/test/suite/codes_tests/test_bse.py b/src/amuse/test/suite/codes_tests/test_bse.py index b5c6f28c5e..5d266a9a05 100644 --- a/src/amuse/test/suite/codes_tests/test_bse.py +++ b/src/amuse/test/suite/codes_tests/test_bse.py @@ -9,7 +9,6 @@ class TestBseInterface(TestWithMPI): - class state(object): def __init__(self): self.type1 = 0.0 @@ -61,26 +60,71 @@ def test1(self): pts2 = 0.01 pts3 = 0.02 sigma = 190.0 - beta = 1.0/8.0 + beta = 1.0 / 8.0 xi = 1.0 - acc2 = 3.0/2.0 + acc2 = 3.0 / 2.0 epsnov = 0.001 eddfac = 1.0 gamma = -1.0 - status = instance.initialize(metallicity, - neta, bwind, hewind, alpha1, CElambda, - ceflag, tflag, ifflag, wdflag, bhflag, - nsflag, mxns, idum, pts1, pts2, pts3, - sigma, beta, xi, acc2, epsnov, eddfac, gamma) + status = instance.initialize( + metallicity, + neta, + bwind, + hewind, + alpha1, + CElambda, + ceflag, + tflag, + ifflag, + wdflag, + bhflag, + nsflag, + mxns, + idum, + pts1, + pts2, + pts3, + sigma, + beta, + xi, + acc2, + epsnov, + eddfac, + gamma, + ) self.assertEqual(status, 0) instance.stop() def test2(self): print("Test basic operations (legacy functions evolve & get_time_step)...") instance = BseInterface() - status = instance.initialize(0.02, 0.5, 0.0, 0.5, 1.0, 0.5, 0, 1, 0, 1, 0, 1, 3.0, - 29769, 0.05, 0.01, 0.02, 190.0, 1.0/8.0, 1.0, 3.0/2.0, 0.001, 1.0, -1.0) + status = instance.initialize( + 0.02, + 0.5, + 0.0, + 0.5, + 1.0, + 0.5, + 0, + 1, + 0, + 1, + 0, + 1, + 3.0, + 29769, + 0.05, + 0.01, + 0.02, + 190.0, + 1.0 / 8.0, + 1.0, + 3.0 / 2.0, + 0.001, + 1.0, + -1.0, + ) new_state = self.state() new_state.mass1 = 3.0 @@ -94,124 +138,309 @@ def test2(self): new_state.eccentricity = 0.5 result = instance.evolve_binary( - new_state.type1, new_state.type2, new_state.initial_mass1, new_state.initial_mass2, - new_state.mass1, new_state.mass2, new_state.radius1, new_state.radius2, - new_state.luminosity1, new_state.luminosity2, new_state.core_mass1, - new_state.core_mass2, new_state.core_radius1, new_state.core_radius2, - new_state.envelope_mass1, new_state.envelope_mass2, new_state.envelope_radius1, - new_state.envelope_radius2, new_state.spin1, new_state.spin2, new_state.epoch1, - new_state.epoch2, new_state.t_ms1, new_state.t_ms2, new_state.bse_age, - new_state.orbital_period, new_state.eccentricity, new_state.end_time + new_state.type1, + new_state.type2, + new_state.initial_mass1, + new_state.initial_mass2, + new_state.mass1, + new_state.mass2, + new_state.radius1, + new_state.radius2, + new_state.luminosity1, + new_state.luminosity2, + new_state.core_mass1, + new_state.core_mass2, + new_state.core_radius1, + new_state.core_radius2, + new_state.envelope_mass1, + new_state.envelope_mass2, + new_state.envelope_radius1, + new_state.envelope_radius2, + new_state.spin1, + new_state.spin2, + new_state.epoch1, + new_state.epoch2, + new_state.t_ms1, + new_state.t_ms2, + new_state.bse_age, + new_state.orbital_period, + new_state.eccentricity, + new_state.end_time, ) updated_state = self.state() - (updated_state.type1, updated_state.type2, updated_state.initial_mass1, updated_state.initial_mass2, - updated_state.mass1, updated_state.mass2, updated_state.radius1, updated_state.radius2, - updated_state.luminosity1, updated_state.luminosity2, updated_state.core_mass1, - updated_state.core_mass2, updated_state.core_radius1, updated_state.core_radius2, - updated_state.envelope_mass1, updated_state.envelope_mass2, updated_state.envelope_radius1, - updated_state.envelope_radius2, updated_state.spin1, updated_state.spin2, - updated_state.epoch1, updated_state.epoch2, updated_state.t_ms1, updated_state.t_ms2, - updated_state.bse_age, updated_state.orbital_period, - updated_state.eccentricity, updated_state.end_time) = result + ( + updated_state.type1, + updated_state.type2, + updated_state.initial_mass1, + updated_state.initial_mass2, + updated_state.mass1, + updated_state.mass2, + updated_state.radius1, + updated_state.radius2, + updated_state.luminosity1, + updated_state.luminosity2, + updated_state.core_mass1, + updated_state.core_mass2, + updated_state.core_radius1, + updated_state.core_radius2, + updated_state.envelope_mass1, + updated_state.envelope_mass2, + updated_state.envelope_radius1, + updated_state.envelope_radius2, + updated_state.spin1, + updated_state.spin2, + updated_state.epoch1, + updated_state.epoch2, + updated_state.t_ms1, + updated_state.t_ms2, + updated_state.bse_age, + updated_state.orbital_period, + updated_state.eccentricity, + updated_state.end_time, + ) = result expected = { - 'radius2': '0x1.c6c8a1c793bcep-1', - 'luminosity2': '0x1.653b1b2d0333bp-1', - 'core_mass2': '0x0.0p+0', - 'bse_age': '0x1.0c6f7a0b5ed8dp-20', - 'end_time': '0x1.0c6f7a0b5ed8dp-20', - 'envelope_mass2': '0x1.0d6fc100ab510p-5', - 'mass2': '0x1.0000000000000p+0', - 'initial_mass2': '0x1.0000000000000p+0', - 'envelope_radius2': '0x1.db27631ba0e5ap-3', - 'core_radius2': '0x0.0p+0', - 'epoch2': '0x0.0p+0', - 't_ms2': '0x1.57d90abe54643p+13', - 'spin2': '0x1.07413b0522aebp+10', + "radius2": "0x1.c6c8a1c793bcep-1", + "luminosity2": "0x1.653b1b2d0333bp-1", + "core_mass2": "0x0.0p+0", + "bse_age": "0x1.0c6f7a0b5ed8dp-20", + "end_time": "0x1.0c6f7a0b5ed8dp-20", + "envelope_mass2": "0x1.0d6fc100ab510p-5", + "mass2": "0x1.0000000000000p+0", + "initial_mass2": "0x1.0000000000000p+0", + "envelope_radius2": "0x1.db27631ba0e5ap-3", + "core_radius2": "0x0.0p+0", + "epoch2": "0x0.0p+0", + "t_ms2": "0x1.57d90abe54643p+13", + "spin2": "0x1.07413b0522aebp+10", } for x in expected: # print "'%s' : '%s'," % (x, getattr(updated_state, x).hex()) - self.assertAlmostRelativeEqual(float.fromhex(expected[x]), getattr(updated_state, x)) + self.assertAlmostRelativeEqual( + float.fromhex(expected[x]), getattr(updated_state, x) + ) self.assertEqual(updated_state.end_time, 1e-06) - dt = instance.get_time_step(updated_state.type1, updated_state.type2, - updated_state.initial_mass1, updated_state.initial_mass2, updated_state.mass1, - updated_state.mass2, updated_state.t_ms1, updated_state.t_ms2, - updated_state.epoch1, updated_state.epoch2, updated_state.bse_age) + dt = instance.get_time_step( + updated_state.type1, + updated_state.type2, + updated_state.initial_mass1, + updated_state.initial_mass2, + updated_state.mass1, + updated_state.mass2, + updated_state.t_ms1, + updated_state.t_ms2, + updated_state.epoch1, + updated_state.epoch2, + updated_state.bse_age, + ) self.assertAlmostEqual(dt, 18.8768, 3) instance.stop() def test3(self): print("Test whether the interface can handle arrays...") instance = BseInterface() - status = instance.initialize(0.02, 0.5, 0.0, 0.5, 1.0, 0.5, 0, 1, 0, 1, 0, 1, 3.0, - 29769, 0.05, 0.01, 0.02, 190.0, 1.0/8.0, 1.0, 3.0/2.0, 0.001, 1.0, -1.0) + status = instance.initialize( + 0.02, + 0.5, + 0.0, + 0.5, + 1.0, + 0.5, + 0, + 1, + 0, + 1, + 0, + 1, + 3.0, + 29769, + 0.05, + 0.01, + 0.02, + 190.0, + 1.0 / 8.0, + 1.0, + 3.0 / 2.0, + 0.001, + 1.0, + -1.0, + ) masses1 = [10.0, 5.0, 4.0] masses2 = [1.0, 1.0, 1.0] types1 = types2 = [1, 1, 1] orbital_periods = [100.0, 200.0, 300.0] eccentricities = [0.5, 0.6, 0.7] - radii1 = luminosity1 = core_mass1 = core_radius1 = envelope_mass1 =\ - envelope_radius1 = spin1 = epoch1 = t_ms1 = [0.0, 0.0, 0.0] - radii2 = luminosity2 = core_mass2 = core_radius2 = envelope_mass2 =\ - envelope_radius2 = spin2 = epoch2 = t_ms2 = [0.0, 0.0, 0.0] + radii1 = ( + luminosity1 + ) = ( + core_mass1 + ) = ( + core_radius1 + ) = envelope_mass1 = envelope_radius1 = spin1 = epoch1 = t_ms1 = [0.0, 0.0, 0.0] + radii2 = ( + luminosity2 + ) = ( + core_mass2 + ) = ( + core_radius2 + ) = envelope_mass2 = envelope_radius2 = spin2 = epoch2 = t_ms2 = [0.0, 0.0, 0.0] init_mass1 = masses1 init_mass2 = masses2 bse_age = [0.0, 0.0, 0.0] end_time = [10.0, 10.0, 10.0] result = instance.evolve_binary( - types1, types2, init_mass1, init_mass2, - masses1, masses2, radii1, radii2, - luminosity1, luminosity2, core_mass1, core_mass2, - core_radius1, core_radius2, envelope_mass1, envelope_mass2, - envelope_radius1, envelope_radius2, spin1, spin2, - epoch1, epoch2, t_ms1, t_ms2, - bse_age, orbital_periods, eccentricities, end_time - ) - self.assertAlmostEqual(result['mass1'][0], 9.977, 2) - self.assertAlmostEqual(result['mass1'][1], 5.0, 2) - self.assertAlmostEqual(result['mass1'][2], 4.0, 2) + types1, + types2, + init_mass1, + init_mass2, + masses1, + masses2, + radii1, + radii2, + luminosity1, + luminosity2, + core_mass1, + core_mass2, + core_radius1, + core_radius2, + envelope_mass1, + envelope_mass2, + envelope_radius1, + envelope_radius2, + spin1, + spin2, + epoch1, + epoch2, + t_ms1, + t_ms2, + bse_age, + orbital_periods, + eccentricities, + end_time, + ) + self.assertAlmostEqual(result["mass1"][0], 9.977, 2) + self.assertAlmostEqual(result["mass1"][1], 5.0, 2) + self.assertAlmostEqual(result["mass1"][2], 4.0, 2) instance.stop() def test4(self): print("Test large number of particles...") number_of_particles = 2000 instance = BseInterface() - status = instance.initialize(0.02, 0.5, 0.0, 0.5, 1.0, 0.5, 0, 1, 0, 1, 0, 1, 3.0, - 29769, 0.05, 0.01, 0.02, 190.0, 1.0/8.0, 1.0, 3.0/2.0, 0.001, 1.0, -1.0) - masses1 = [1.0 + ((x / 1.0*number_of_particles) * 10.0) for x in range(1, number_of_particles+1)] - masses2 = [2.0 + ((x / 1.0*number_of_particles) * 5.0) for x in range(1, number_of_particles+1)] - orbital_periods = [100.0 + ((x / 1.0*number_of_particles) * 900.0) for x in range(1, number_of_particles+1)] - eccentricities = [0.5 + ((x / 1.0*number_of_particles) * 0.4) for x in range(1, number_of_particles+1)] - - types1 = types2 = [1 for x in range(1, number_of_particles+1)] - radii1 = luminosity1 = core_mass1 = core_radius1 = envelope_mass1 =\ - envelope_radius1 = spin1 = epoch1 = t_ms1 =\ - radii2 = luminosity2 = core_mass2 = core_radius2 = envelope_mass2 =\ - envelope_radius2 = spin2 = epoch2 = t_ms2 =\ - bse_age = [0.0 for x in range(1, number_of_particles+1)] - end_time = [1.0 for x in range(1, number_of_particles+1)] + status = instance.initialize( + 0.02, + 0.5, + 0.0, + 0.5, + 1.0, + 0.5, + 0, + 1, + 0, + 1, + 0, + 1, + 3.0, + 29769, + 0.05, + 0.01, + 0.02, + 190.0, + 1.0 / 8.0, + 1.0, + 3.0 / 2.0, + 0.001, + 1.0, + -1.0, + ) + masses1 = [ + 1.0 + ((x / 1.0 * number_of_particles) * 10.0) + for x in range(1, number_of_particles + 1) + ] + masses2 = [ + 2.0 + ((x / 1.0 * number_of_particles) * 5.0) + for x in range(1, number_of_particles + 1) + ] + orbital_periods = [ + 100.0 + ((x / 1.0 * number_of_particles) * 900.0) + for x in range(1, number_of_particles + 1) + ] + eccentricities = [ + 0.5 + ((x / 1.0 * number_of_particles) * 0.4) + for x in range(1, number_of_particles + 1) + ] + + types1 = types2 = [1 for x in range(1, number_of_particles + 1)] + radii1 = ( + luminosity1 + ) = ( + core_mass1 + ) = ( + core_radius1 + ) = ( + envelope_mass1 + ) = ( + envelope_radius1 + ) = ( + spin1 + ) = ( + epoch1 + ) = ( + t_ms1 + ) = ( + radii2 + ) = ( + luminosity2 + ) = ( + core_mass2 + ) = ( + core_radius2 + ) = envelope_mass2 = envelope_radius2 = spin2 = epoch2 = t_ms2 = bse_age = [ + 0.0 for x in range(1, number_of_particles + 1) + ] + end_time = [1.0 for x in range(1, number_of_particles + 1)] init_mass1 = masses1 init_mass2 = masses2 result = instance.evolve_binary( - types1, types2, init_mass1, init_mass2, - masses1, masses2, radii1, radii2, - luminosity1, luminosity2, core_mass1, core_mass2, - core_radius1, core_radius2, envelope_mass1, envelope_mass2, - envelope_radius1, envelope_radius2, spin1, spin2, - epoch1, epoch2, t_ms1, t_ms2, - bse_age, orbital_periods, eccentricities, end_time - ) - self.assertEqual(len(result['mass1']), number_of_particles) + types1, + types2, + init_mass1, + init_mass2, + masses1, + masses2, + radii1, + radii2, + luminosity1, + luminosity2, + core_mass1, + core_mass2, + core_radius1, + core_radius2, + envelope_mass1, + envelope_mass2, + envelope_radius1, + envelope_radius2, + spin1, + spin2, + epoch1, + epoch2, + t_ms1, + t_ms2, + bse_age, + orbital_periods, + eccentricities, + end_time, + ) + self.assertEqual(len(result["mass1"]), number_of_particles) instance.stop() class TestBse(TestWithMPI): - def test1(self): print("Testing evolution of a close binary system...") instance = Bse() @@ -225,7 +454,9 @@ def test1(self): stars[1].mass = 0.3 | units.MSun orbital_period = 200.0 | units.day - semi_major_axis = instance.orbital_period_to_semi_major_axis(orbital_period, stars[0].mass, stars[1].mass) + semi_major_axis = instance.orbital_period_to_semi_major_axis( + orbital_period, stars[0].mass, stars[1].mass + ) instance.particles.add_particles(stars) @@ -252,12 +483,16 @@ def test1(self): while current_time < (480 | units.Myr): instance.update_time_steps() # The next line appears a bit weird, but saves time for this simple test. - current_time = current_time + max(5.0*instance.binaries[0].time_step, 0.3 | units.Myr) + current_time = current_time + max( + 5.0 * instance.binaries[0].time_step, 0.3 | units.Myr + ) instance.evolve_model(current_time) from_bse_to_model.copy() from_bse_to_model_binaries.copy() if not binary.child1.stellar_type == previous_type: - results.append((binary.age, binary.child1.mass, binary.child1.stellar_type)) + results.append( + (binary.age, binary.child1.mass, binary.child1.stellar_type) + ) previous_type = binary.child1.stellar_type self.assertEqual(len(results), 6) @@ -283,7 +518,9 @@ def test1(self): 332.2786 | units.Myr, ) for result, expected in zip(results, times): - self.assertAlmostEqual(result[0].value_in(units.Myr), expected.value_in(units.Myr), 0) + self.assertAlmostEqual( + result[0].value_in(units.Myr), expected.value_in(units.Myr), 0 + ) masses = ( 3.000 | units.MSun, @@ -294,7 +531,9 @@ def test1(self): 0.707 | units.MSun, ) for result, expected in zip(results, masses): - self.assertAlmostEqual(result[1].value_in(units.MSun), expected.value_in(units.MSun), 2) + self.assertAlmostEqual( + result[1].value_in(units.MSun), expected.value_in(units.MSun), 2 + ) instance.stop() @@ -310,7 +549,9 @@ def test2(self): stars[0].mass = 3.0 | units.MSun stars[1].mass = 0.3 | units.MSun orbital_period = 2.0e5 | units.day - semi_major_axis = instance.orbital_period_to_semi_major_axis(orbital_period, stars[0].mass, stars[1].mass) + semi_major_axis = instance.orbital_period_to_semi_major_axis( + orbital_period, stars[0].mass, stars[1].mass + ) instance.particles.add_particles(stars) @@ -337,12 +578,16 @@ def test2(self): while current_time < (335 | units.Myr): instance.update_time_steps() # The next line appears a bit weird, but saves time for this simple test. - current_time = current_time + max(2.0*instance.binaries[0].time_step, 0.04 | units.Myr) + current_time = current_time + max( + 2.0 * instance.binaries[0].time_step, 0.04 | units.Myr + ) instance.evolve_model(current_time) from_bse_to_model.copy() from_bse_to_model_binaries.copy() if not binary.child1.stellar_type == previous_type: - results.append((binary.age, binary.child1.mass, binary.child1.stellar_type)) + results.append( + (binary.age, binary.child1.mass, binary.child1.stellar_type) + ) previous_type = binary.child1.stellar_type print(results) self.assertEqual(len(results), 6) @@ -353,10 +598,12 @@ def test2(self): 287.7848 | units.Myr, 331.1454 | units.Myr, 332.7407 | units.Myr, - 333.4146 | units.Myr + 333.4146 | units.Myr, ) for result, expected in zip(results, times): - self.assertAlmostEqual(result[0].value_in(units.Myr), expected.value_in(units.Myr), 0) + self.assertAlmostEqual( + result[0].value_in(units.Myr), expected.value_in(units.Myr), 0 + ) masses = ( 3.000 | units.MSun, @@ -364,10 +611,12 @@ def test2(self): 2.999 | units.MSun, 2.956 | units.MSun, 2.919 | units.MSun, - 0.928 | units.MSun + 0.928 | units.MSun, ) for result, expected in zip(results, masses): - self.assertAlmostEqual(result[1].value_in(units.MSun), expected.value_in(units.MSun), 2) + self.assertAlmostEqual( + result[1].value_in(units.MSun), expected.value_in(units.MSun), 2 + ) types = ( "Hertzsprung Gap", @@ -395,7 +644,9 @@ def test3(self): stars[1].mass = 4.387 | units.MSun orbital_period = 1964.18453 | units.day - semi_major_axis = instance.orbital_period_to_semi_major_axis(orbital_period, stars[0].mass, stars[1].mass) + semi_major_axis = instance.orbital_period_to_semi_major_axis( + orbital_period, stars[0].mass, stars[1].mass + ) instance.particles.add_particles(stars) binaries = Particles(1) @@ -422,16 +673,28 @@ def test3(self): while current_time < (170 | units.Myr): instance.update_time_steps() # The next line appears a bit weird, but saves time for this simple test. - current_time = current_time + max(2.0*instance.binaries[0].time_step, 0.04 | units.Myr) + current_time = current_time + max( + 2.0 * instance.binaries[0].time_step, 0.04 | units.Myr + ) instance.evolve_model(current_time) from_bse_to_model.copy() from_bse_to_model_binaries.copy() - if not (binary.child1.stellar_type == previous_type1 and binary.child2.stellar_type == previous_type2): - results.append((binary.age, str(binary.child1.stellar_type)+" and "+str(binary.child2.stellar_type))) + if not ( + binary.child1.stellar_type == previous_type1 + and binary.child2.stellar_type == previous_type2 + ): + results.append( + ( + binary.age, + str(binary.child1.stellar_type) + + " and " + + str(binary.child2.stellar_type), + ) + ) previous_type1 = binary.child1.stellar_type previous_type2 = binary.child2.stellar_type - print('\n'.join(map(str, results))) + print("\n".join(map(str, results))) self.assertEqual(len(results), 13) times = ( 38.9708 | units.Myr, @@ -446,10 +709,12 @@ def test3(self): 166.1043 | units.Myr, 166.5795 | units.Myr, 166.9627 | units.Myr, - 166.9863 | units.Myr + 166.9863 | units.Myr, ) for result, expected in zip(results, times): - self.assertAlmostEqual(result[0].value_in(units.Myr), expected.value_in(units.Myr), 0) + self.assertAlmostEqual( + result[0].value_in(units.Myr), expected.value_in(units.Myr), 0 + ) types = ( "Hertzsprung Gap and Main Sequence star", @@ -492,7 +757,9 @@ def test4(self): binary = binaries[0] orbital_period = 1964.18453 | units.day - semi_major_axis = instance.orbital_period_to_semi_major_axis(orbital_period, stars[0].mass, stars[1].mass) + semi_major_axis = instance.orbital_period_to_semi_major_axis( + orbital_period, stars[0].mass, stars[1].mass + ) binary.semi_major_axis = semi_major_axis binary.eccentricity = 0.0 binary.child1 = stars[0] @@ -534,7 +801,9 @@ def test5(self): binary = binaries[0] orbital_period = 200.0 | units.day - semi_major_axis = instance.orbital_period_to_semi_major_axis(orbital_period, stars[0].mass, stars[1].mass) + semi_major_axis = instance.orbital_period_to_semi_major_axis( + orbital_period, stars[0].mass, stars[1].mass + ) binary.semi_major_axis = semi_major_axis binary.eccentricity = 0.99 binary.child1 = stars[0] @@ -594,12 +863,10 @@ def test7(self): binaries.eccentricity = 0.0 for i in range(3): binaries[i].child1 = stars[i] - binaries[i].child2 = stars[i+3] + binaries[i].child2 = stars[i + 3] orbital_period = 200.0 | units.day semi_major_axis = instance.orbital_period_to_semi_major_axis( - orbital_period, - binaries.child1.as_set().mass, - binaries.child2.as_set().mass + orbital_period, binaries.child1.as_set().mass, binaries.child2.as_set().mass ) binaries.semi_major_axis = semi_major_axis @@ -607,24 +874,46 @@ def test7(self): instance.binaries.add_particles(binaries) self.assertAlmostEqual(instance.binaries.age, [0.0, 0.0, 0.0] | units.yr) - self.assertAlmostEqual(instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) - print("evolve_model without arguments: use shared timestep = min(particles.time_step)") + print( + "evolve_model without arguments: use shared timestep = min(particles.time_step)" + ) instance.evolve_model() - self.assertAlmostEqual(instance.binaries.age, [18.8768, 18.8768, 18.8768] | units.Myr, 3) - self.assertAlmostEqual(instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.binaries.age, [18.8768, 18.8768, 18.8768] | units.Myr, 3 + ) + self.assertAlmostEqual( + instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) self.assertAlmostEqual(instance.model_time, 18.8768 | units.Myr, 3) - print("evolve_model with end_time: take timesteps, until end_time is reached exactly") + print( + "evolve_model with end_time: take timesteps, until end_time is reached exactly" + ) instance.evolve_model(100 | units.Myr) - self.assertAlmostEqual(instance.binaries.age, [100.0, 100.0, 100.0] | units.Myr, 3) - self.assertAlmostEqual(instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.binaries.age, [100.0, 100.0, 100.0] | units.Myr, 3 + ) + self.assertAlmostEqual( + instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) self.assertAlmostEqual(instance.model_time, 100.0 | units.Myr, 3) - print("evolve_model with keep_synchronous: use non-shared timestep, particle ages will typically diverge") + print( + "evolve_model with keep_synchronous: use non-shared timestep, particle ages will typically diverge" + ) instance.evolve_model(keep_synchronous=False) - self.assertAlmostEqual(instance.binaries.age, (100 | units.Myr) + ([550.1565, 58.2081, 18.8768] | units.Myr), 3) - self.assertAlmostEqual(instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.binaries.age, + (100 | units.Myr) + ([550.1565, 58.2081, 18.8768] | units.Myr), + 3, + ) + self.assertAlmostEqual( + instance.binaries.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) self.assertAlmostEqual(instance.model_time, 100.0 | units.Myr, 3) # Unchanged! instance.stop() @@ -641,12 +930,10 @@ def test8(self): binaries.eccentricity = 0.0 for i in range(3): binaries[i].child1 = stars[i] - binaries[i].child2 = stars[i+3] + binaries[i].child2 = stars[i + 3] orbital_period = 200.0 | units.day semi_major_axis = instance.orbital_period_to_semi_major_axis( - orbital_period, - binaries.child1.as_set().mass, - binaries.child2.as_set().mass + orbital_period, binaries.child1.as_set().mass, binaries.child2.as_set().mass ) binaries.semi_major_axis = semi_major_axis @@ -669,9 +956,13 @@ def test8(self): self.assertEqual(len(instance.binaries), 3) # it's back... self.assertAlmostEqual(instance.binaries[0].age, 2.0 | units.Myr) self.assertAlmostEqual(instance.binaries[1].age, 0.0 | units.Myr) - self.assertAlmostEqual(instance.binaries[2].age, 0.0 | units.Myr) # ... and rejuvenated. + self.assertAlmostEqual( + instance.binaries[2].age, 0.0 | units.Myr + ) # ... and rejuvenated. - instance.evolve_model(3.0 | units.Myr) # The young stars keep their age offset from the old star + instance.evolve_model( + 3.0 | units.Myr + ) # The young stars keep their age offset from the old star self.assertAlmostEqual(instance.binaries.age, [3.0, 1.0, 1.0] | units.Myr) instance.evolve_model(4.0 | units.Myr) self.assertAlmostEqual(instance.binaries.age, [4.0, 2.0, 2.0] | units.Myr) @@ -686,32 +977,37 @@ def test9(self): binaries = Particles(1) orbital_period = 200.0 | units.day - semi_major_axis = instance.orbital_period_to_semi_major_axis(orbital_period, stars[0].mass, stars[1].mass) + semi_major_axis = instance.orbital_period_to_semi_major_axis( + orbital_period, stars[0].mass, stars[1].mass + ) binaries.semi_major_axis = semi_major_axis binaries.eccentricity = 0.0 binaries[0].child1 = stars[0] binaries[0].child2 = stars[1] - print("First do everything manually:", end=' ') - self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') + print("First do everything manually:", end=" ") + self.assertEqual(instance.get_name_of_current_state(), "UNINITIALIZED") instance.initialize_code() - self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "INITIALIZED") instance.commit_parameters() - self.assertEqual(instance.get_name_of_current_state(), 'RUN') + self.assertEqual(instance.get_name_of_current_state(), "RUN") instance.cleanup_code() - self.assertEqual(instance.get_name_of_current_state(), 'END') + self.assertEqual(instance.get_name_of_current_state(), "END") instance.stop() print("ok") - print("initialize_code(), commit_parameters(), " - "and cleanup_code() should be called automatically:", end=' ') + print( + "initialize_code(), commit_parameters(), " + "and cleanup_code() should be called automatically:", + end=" ", + ) instance = Bse() - self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "UNINITIALIZED") instance.parameters.reimers_mass_loss_coefficient = 0.5 - self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "INITIALIZED") instance.particles.add_particles(stars) instance.binaries.add_particles(binaries) - self.assertEqual(instance.get_name_of_current_state(), 'RUN') + self.assertEqual(instance.get_name_of_current_state(), "RUN") instance.stop() - self.assertEqual(instance.get_name_of_current_state(), 'STOPPED') + self.assertEqual(instance.get_name_of_current_state(), "STOPPED") print("ok") diff --git a/src/amuse/test/suite/codes_tests/test_sse.py b/src/amuse/test/suite/codes_tests/test_sse.py index 0fbde2b209..2701f3c48a 100644 --- a/src/amuse/test/suite/codes_tests/test_sse.py +++ b/src/amuse/test/suite/codes_tests/test_sse.py @@ -15,7 +15,6 @@ class TestMPIInterface(TestWithMPI): - class state(object): def __init__(self): self.stellar_type = 0.0 @@ -50,10 +49,21 @@ def initialize_module_with_default_parameters(self, sse): pts2 = 0.01 pts3 = 0.02 - status = sse.initialize(metallicity, - neta, bwind, hewind, sigma, - ifflag, wdflag, bhflag, nsflag, mxns, - pts1, pts2, pts3) + status = sse.initialize( + metallicity, + neta, + bwind, + hewind, + sigma, + ifflag, + wdflag, + bhflag, + nsflag, + mxns, + pts1, + pts2, + pts3, + ) def test1(self): sse = SseInterface() @@ -75,10 +85,21 @@ def test1(self): pts2 = 0.01 pts3 = 0.02 - status = sse.initialize(metallicity, - neta, bwind, hewind, sigma, - ifflag, wdflag, bhflag, nsflag, mxns, - pts1, pts2, pts3) + status = sse.initialize( + metallicity, + neta, + bwind, + hewind, + sigma, + ifflag, + wdflag, + bhflag, + nsflag, + mxns, + pts1, + pts2, + pts3, + ) self.assertEqual(status, 0) @@ -104,10 +125,21 @@ def test2(self): pts2 = 0.01 pts3 = 0.02 - status = sse.initialize(metallicity, - neta, bwind, hewind, sigma, - ifflag, wdflag, bhflag, nsflag, mxns, - pts1, pts2, pts3) + status = sse.initialize( + metallicity, + neta, + bwind, + hewind, + sigma, + ifflag, + wdflag, + bhflag, + nsflag, + mxns, + pts1, + pts2, + pts3, + ) self.assertEqual(status, 0) new_state = self.state() new_state.mass = 1.0 @@ -116,45 +148,84 @@ def test2(self): new_state.age = 1e-06 result = sse.evolve_star( new_state.stellar_type, - new_state.zams_mass, new_state.mass, new_state.radius, - new_state.luminosity, new_state.core_mass, new_state.core_radius, - new_state.envelope_mass, new_state.envelope_radius, new_state.spin, - new_state.epoch, new_state.t_ms, new_state.sse_age, new_state.age + new_state.zams_mass, + new_state.mass, + new_state.radius, + new_state.luminosity, + new_state.core_mass, + new_state.core_radius, + new_state.envelope_mass, + new_state.envelope_radius, + new_state.spin, + new_state.epoch, + new_state.t_ms, + new_state.sse_age, + new_state.age, ) updated_state = self.state() - (updated_state.stellar_type, updated_state.zams_mass, updated_state.mass, updated_state.radius, - updated_state.luminosity, updated_state.core_mass, updated_state.core_radius, - updated_state.envelope_mass, updated_state.envelope_radius, updated_state.spin, - updated_state.epoch, updated_state.t_ms, updated_state.sse_age, updated_state.age) = result - attributes = ('stellar_type', 'zams_mass', 'mass', 'radius', 'luminosity', 'core_mass', 'core_radius', - 'envelope_mass', 'envelope_radius', 'spin', 'epoch', 't_ms', 'sse_age', 'age') + ( + updated_state.stellar_type, + updated_state.zams_mass, + updated_state.mass, + updated_state.radius, + updated_state.luminosity, + updated_state.core_mass, + updated_state.core_radius, + updated_state.envelope_mass, + updated_state.envelope_radius, + updated_state.spin, + updated_state.epoch, + updated_state.t_ms, + updated_state.sse_age, + updated_state.age, + ) = result + attributes = ( + "stellar_type", + "zams_mass", + "mass", + "radius", + "luminosity", + "core_mass", + "core_radius", + "envelope_mass", + "envelope_radius", + "spin", + "epoch", + "t_ms", + "sse_age", + "age", + ) expected = { - 'zams_mass': '0x1.0000000000000p+0', - 'mass': '0x1.0000000000000p+0', - 'radius': '0x1.c6c8a1c793bcep-1', - 'luminosity': '0x1.653b1b2d0333bp-1', - 'core_mass': '0x0.0p+0', - 'core_radius': '0x0.0p+0', - 'envelope_mass': '0x1.0d6fc100ab510p-5', - 'envelope_radius': '0x1.db27631ba0e5ap-3', - 'spin': '0x1.07413b0522d63p+10', - 'epoch': '0x0.0p+0', - 't_ms': '0x1.57d90abe54643p+13', - 'sse_age': '0x1.0c6f7a0b5ed8dp-20', - 'age': '0x1.0c6f7a0b5ed8dp-20', + "zams_mass": "0x1.0000000000000p+0", + "mass": "0x1.0000000000000p+0", + "radius": "0x1.c6c8a1c793bcep-1", + "luminosity": "0x1.653b1b2d0333bp-1", + "core_mass": "0x0.0p+0", + "core_radius": "0x0.0p+0", + "envelope_mass": "0x1.0d6fc100ab510p-5", + "envelope_radius": "0x1.db27631ba0e5ap-3", + "spin": "0x1.07413b0522d63p+10", + "epoch": "0x0.0p+0", + "t_ms": "0x1.57d90abe54643p+13", + "sse_age": "0x1.0c6f7a0b5ed8dp-20", + "age": "0x1.0c6f7a0b5ed8dp-20", } for x in expected: - self.assertAlmostRelativeEqual(float.fromhex(expected[x]), getattr(updated_state, x)) + self.assertAlmostRelativeEqual( + float.fromhex(expected[x]), getattr(updated_state, x) + ) self.assertEqual(updated_state.age, 1e-06) - dt = sse.get_time_step(updated_state.stellar_type, + dt = sse.get_time_step( + updated_state.stellar_type, updated_state.zams_mass, updated_state.age, updated_state.mass, updated_state.t_ms, - updated_state.epoch) + updated_state.epoch, + ) self.assertAlmostEqual(dt, 550.1565, 2) sse.stop() @@ -164,8 +235,13 @@ def test3(self): types = [1, 1, 1] masses = [10, 5, 4] radii = [5.0, 2.0, 1.0] - luminosity = core_mass = core_radius = envelope_mass =\ - envelope_radius = spin = epoch = t_ms = [0.0, 0.0, 0.0] + luminosity = ( + core_mass + ) = core_radius = envelope_mass = envelope_radius = spin = epoch = t_ms = [ + 0.0, + 0.0, + 0.0, + ] sse_age = age = [1e-6, 1e-06, 1e-6] result = sse.evolve_star( types, @@ -181,11 +257,11 @@ def test3(self): epoch, t_ms, sse_age, - age + age, ) - self.assertEqual(result['mass'][0], 10) - self.assertEqual(result['mass'][1], 5) - self.assertAlmostEqual(result['mass'][2], 4.0, 2) + self.assertEqual(result["mass"][0], 10) + self.assertEqual(result["mass"][1], 5) + self.assertAlmostEqual(result["mass"][2], 4.0, 2) sse.stop() def test4(self): @@ -194,9 +270,11 @@ def test4(self): types = [1 for x in range(1, 4000)] masses = [1.0 + ((x / 4000.0) * 10.0) for x in range(1, 4000)] radii = [1.0 for x in range(1, 4000)] - luminosity = core_mass = core_radius = envelope_mass =\ - envelope_radius = spin = epoch =\ - t_ms = [0.0 for x in range(1, 4000)] + luminosity = ( + core_mass + ) = core_radius = envelope_mass = envelope_radius = spin = epoch = t_ms = [ + 0.0 for x in range(1, 4000) + ] sse_age = age = [1e-06 for x in range(1, 4000)] result = sse.evolve_star( @@ -213,14 +291,13 @@ def test4(self): epoch, t_ms, sse_age, - age + age, ) - self.assertEqual(len(result['mass']), 3999) + self.assertEqual(len(result["mass"]), 3999) sse.stop() class TestSse(TestWithMPI): - def test1(self): sse = Sse() sse.commit_parameters() @@ -259,10 +336,12 @@ def test1(self): 104.7 | units.Myr, 120.1 | units.Myr, 120.9 | units.Myr, - 121.5 | units.Myr + 121.5 | units.Myr, ) for result, expected in zip(results, times): - self.assertAlmostEqual(result[0].value_in(units.Myr), expected.value_in(units.Myr), 1) + self.assertAlmostEqual( + result[0].value_in(units.Myr), expected.value_in(units.Myr), 1 + ) masses = ( 5.000 | units.MSun, @@ -270,10 +349,12 @@ def test1(self): 4.998 | units.MSun, 4.932 | units.MSun, 4.895 | units.MSun, - 0.997 | units.MSun + 0.997 | units.MSun, ) for result, expected in zip(results, masses): - self.assertAlmostEqual(result[1].value_in(units.MSun), expected.value_in(units.MSun), 3) + self.assertAlmostEqual( + result[1].value_in(units.MSun), expected.value_in(units.MSun), 3 + ) types = ( "Hertzsprung Gap", @@ -302,7 +383,9 @@ def test2(self): sse.evolve_model(120.1 | units.Myr) self.assertAlmostEqual(sse.particles[0].mass.value_in(units.MSun), 4.932, 3) - self.assertAlmostEqual(sse.particles[0].temperature.value_in(units.K), 4221., 0) + self.assertAlmostEqual( + sse.particles[0].temperature.value_in(units.K), 4221.0, 0 + ) sse.stop() @@ -360,13 +443,13 @@ def test5(self): def test6(self): print("Test whether a set of stars evolves synchronously...") -# Create an array of stars with a range in stellar mass - masses = [.5, 1., 2., 5., 10., 30.] | units.MSun + # Create an array of stars with a range in stellar mass + masses = [0.5, 1.0, 2.0, 5.0, 10.0, 30.0] | units.MSun number_of_stars = len(masses) stars = Particles(number_of_stars) stars.mass = masses -# Initialize stellar evolution code + # Initialize stellar evolution code instance = Sse() instance.commit_parameters() instance.particles.add_particles(stars) @@ -394,7 +477,9 @@ def test6(self): def test7(self): print("Test: evolve particles one at a time.") - print("Used to be problematic, since initial_mass of idle particle is set to zero.") + print( + "Used to be problematic, since initial_mass of idle particle is set to zero." + ) stars = Particles(2) stars.mass = 1.0 | units.MSun for star in stars: @@ -402,7 +487,9 @@ def test7(self): stellar_evolution.commit_parameters() stellar_evolution.particles.add_particles(star.as_set()) stellar_evolution.commit_particles() - from_stellar_evolution_to_model = stellar_evolution.particles.new_channel_to(star.as_set()) + from_stellar_evolution_to_model = ( + stellar_evolution.particles.new_channel_to(star.as_set()) + ) stellar_evolution.evolve_model() from_stellar_evolution_to_model.copy() stellar_evolution.stop() @@ -449,14 +536,16 @@ def test10(self): stars = Particles(10) stars.mass = 1.0 | units.MSun stellar_evolution.particles.add_particles(stars) - self.assertEqual(stellar_evolution.particles._factory_for_new_collection(), Particles) + self.assertEqual( + stellar_evolution.particles._factory_for_new_collection(), Particles + ) filename = os.path.join(get_path_to_results(), "test.h5") if os.path.exists(filename): os.remove(filename) - io.write_set_to_file(stellar_evolution.particles, filename, 'hdf5') - stored_stars = io.read_set_from_file(filename, 'hdf5') + io.write_set_to_file(stellar_evolution.particles, filename, "hdf5") + stored_stars = io.read_set_from_file(filename, "hdf5") self.assertEqual(len(stars), len(stored_stars)) self.assertAlmostRelativeEquals(stars.mass, stored_stars.mass) @@ -470,25 +559,50 @@ def test11(self): instance.particles.add_particles(stars) self.assertEqual(instance.particles.age, [0.0, 0.0, 0.0] | units.yr) - self.assertAlmostEqual(instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) - self.assertAlmostEqual(instance.particles.radius, [0.8882494502, 1.610210385, 1.979134445] | units.RSun) + self.assertAlmostEqual( + instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) + self.assertAlmostEqual( + instance.particles.radius, + [0.8882494502, 1.610210385, 1.979134445] | units.RSun, + ) - print("evolve_model without arguments: use shared timestep = min(particles.time_step)") + print( + "evolve_model without arguments: use shared timestep = min(particles.time_step)" + ) instance.evolve_model() - self.assertAlmostEqual(instance.particles.age, [18.8768, 18.8768, 18.8768] | units.Myr, 3) - self.assertAlmostEqual(instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.particles.age, [18.8768, 18.8768, 18.8768] | units.Myr, 3 + ) + self.assertAlmostEqual( + instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) self.assertAlmostEqual(instance.model_time, 18.8768 | units.Myr, 3) - print("evolve_model with end_time: take timesteps, until end_time is reached exactly") + print( + "evolve_model with end_time: take timesteps, until end_time is reached exactly" + ) instance.evolve_model(100 | units.Myr) - self.assertAlmostEqual(instance.particles.age, [100.0, 100.0, 100.0] | units.Myr, 3) - self.assertAlmostEqual(instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.particles.age, [100.0, 100.0, 100.0] | units.Myr, 3 + ) + self.assertAlmostEqual( + instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) self.assertAlmostEqual(instance.model_time, 100.0 | units.Myr, 3) - print("evolve_model with keep_synchronous: use non-shared timestep, particle ages will typically diverge") + print( + "evolve_model with keep_synchronous: use non-shared timestep, particle ages will typically diverge" + ) instance.evolve_model(keep_synchronous=False) - self.assertAlmostEqual(instance.particles.age, (100 | units.Myr) + ([550.1565, 58.2081, 18.8768] | units.Myr), 3) - self.assertAlmostEqual(instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3) + self.assertAlmostEqual( + instance.particles.age, + (100 | units.Myr) + ([550.1565, 58.2081, 18.8768] | units.Myr), + 3, + ) + self.assertAlmostEqual( + instance.particles.time_step, [550.1565, 58.2081, 18.8768] | units.Myr, 3 + ) self.assertAlmostEqual(instance.model_time, 100.0 | units.Myr, 3) # Unchanged! instance.stop() @@ -517,9 +631,13 @@ def test12(self): self.assertEqual(len(instance.particles), 3) # it's back... self.assertAlmostEqual(instance.particles[0].age, 2.0 | units.Myr) self.assertAlmostEqual(instance.particles[1].age, 0.0 | units.Myr) - self.assertAlmostEqual(instance.particles[2].age, 0.0 | units.Myr) # ... and rejuvenated. + self.assertAlmostEqual( + instance.particles[2].age, 0.0 | units.Myr + ) # ... and rejuvenated. - instance.evolve_model(3.0 | units.Myr) # The young stars keep their age offset from the old star + instance.evolve_model( + 3.0 | units.Myr + ) # The young stars keep their age offset from the old star self.assertAlmostEqual(instance.particles.age, [3.0, 1.0, 1.0] | units.Myr) instance.evolve_model(4.0 | units.Myr) self.assertAlmostEqual(instance.particles.age, [4.0, 2.0, 2.0] | units.Myr) @@ -530,28 +648,31 @@ def test13(self): stars = Particles(1) stars.mass = 1.0 | units.MSun - print("First do everything manually:", end=' ') + print("First do everything manually:", end=" ") instance = Sse() - self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "UNINITIALIZED") instance.initialize_code() - self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "INITIALIZED") instance.commit_parameters() - self.assertEqual(instance.get_name_of_current_state(), 'RUN') + self.assertEqual(instance.get_name_of_current_state(), "RUN") instance.cleanup_code() - self.assertEqual(instance.get_name_of_current_state(), 'END') + self.assertEqual(instance.get_name_of_current_state(), "END") instance.stop() print("ok") - print("initialize_code(), commit_parameters(), " - "and cleanup_code() should be called automatically:", end=' ') + print( + "initialize_code(), commit_parameters(), " + "and cleanup_code() should be called automatically:", + end=" ", + ) instance = Sse() - self.assertEqual(instance.get_name_of_current_state(), 'UNINITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "UNINITIALIZED") instance.parameters.reimers_mass_loss_coefficient = 0.5 - self.assertEqual(instance.get_name_of_current_state(), 'INITIALIZED') + self.assertEqual(instance.get_name_of_current_state(), "INITIALIZED") instance.particles.add_particles(stars) - self.assertEqual(instance.get_name_of_current_state(), 'RUN') + self.assertEqual(instance.get_name_of_current_state(), "RUN") instance.stop() - self.assertEqual(instance.get_name_of_current_state(), 'STOPPED') + self.assertEqual(instance.get_name_of_current_state(), "STOPPED") print("ok") def test14a(self): @@ -570,10 +691,14 @@ def test14a(self): for i in range(1, number_of_steps + 1): se_stars[1].evolve_for(step_size) self.assertAlmostEqual(se_stars.age, [number_of_steps, i] * step_size) - self.assertAlmostRelativeEqual(se_stars[0].age, se_stars[1].age) - self.assertAlmostRelativeEqual(se_stars[0].luminosity, se_stars[1].luminosity, 3) - self.assertAlmostRelativeEqual(se_stars[0].radius, se_stars[1].radius, 3) - self.assertAlmostRelativeEqual(se_stars[0].temperature, se_stars[1].temperature, 3) + self.assertAlmostRelativeEqual(se_stars[0].age, se_stars[1].age) + self.assertAlmostRelativeEqual( + se_stars[0].luminosity, se_stars[1].luminosity, 3 + ) + self.assertAlmostRelativeEqual(se_stars[0].radius, se_stars[1].radius, 3) + self.assertAlmostRelativeEqual( + se_stars[0].temperature, se_stars[1].temperature, 3 + ) instance.stop() def test14b(self): @@ -592,10 +717,14 @@ def test14b(self): for i in range(1, number_of_steps + 1): se_stars[1:].evolve_for(step_size) self.assertAlmostEqual(se_stars.age, [number_of_steps, i] * step_size) - self.assertAlmostRelativeEqual(se_stars[0].age, se_stars[1].age) - self.assertAlmostRelativeEqual(se_stars[0].luminosity, se_stars[1].luminosity, 3) - self.assertAlmostRelativeEqual(se_stars[0].radius, se_stars[1].radius, 3) - self.assertAlmostRelativeEqual(se_stars[0].temperature, se_stars[1].temperature, 3) + self.assertAlmostRelativeEqual(se_stars[0].age, se_stars[1].age) + self.assertAlmostRelativeEqual( + se_stars[0].luminosity, se_stars[1].luminosity, 3 + ) + self.assertAlmostRelativeEqual(se_stars[0].radius, se_stars[1].radius, 3) + self.assertAlmostRelativeEqual( + se_stars[0].temperature, se_stars[1].temperature, 3 + ) instance.stop() def test15(self): @@ -605,16 +734,17 @@ def test15(self): class notsorandom(object): def random(self, N): - return numpy.array(range(N))/(N-1.) + return numpy.array(range(N)) / (N - 1.0) def random_sample(self, N): - return numpy.array(range(N))/(N-1.) + return numpy.array(range(N)) / (N - 1.0) masses = new_salpeter_mass_distribution( number_of_stars, mass_min=0.1 | units.MSun, mass_max=100.0 | units.MSun, - alpha=-1.01, random=notsorandom() + alpha=-1.01, + random=notsorandom(), ) stars = Particles(mass=masses) @@ -635,16 +765,17 @@ def test16(self): class notsorandom(object): def random(self, N): - return numpy.array(range(N))/(N-1.) + return numpy.array(range(N)) / (N - 1.0) def random_sample(self, N): - return numpy.array(range(N))/(N-1.) + return numpy.array(range(N)) / (N - 1.0) masses = new_salpeter_mass_distribution( number_of_stars, mass_min=0.1 | units.MSun, mass_max=100.0 | units.MSun, - alpha=-1.01, random=notsorandom() + alpha=-1.01, + random=notsorandom(), ) stars = Particles(mass=masses) @@ -675,8 +806,10 @@ def test17(self): revived.evolve_for(numpy.pi | units.Myr) for star in instance.particles: star.evolve_for(star.age) - self.assertAlmostEqual(instance.particles.age[:-1], 2*time_steps[[0, 2, 3, 5, 6, 7, 9]]) - self.assertAlmostEqual(instance.particles.age[-1], 2*numpy.pi | units.Myr) + self.assertAlmostEqual( + instance.particles.age[:-1], 2 * time_steps[[0, 2, 3, 5, 6, 7, 9]] + ) + self.assertAlmostEqual(instance.particles.age[-1], 2 * numpy.pi | units.Myr) instance.particles.remove_particles(particles[[2, 5, 6]]) instance.particles.add_particles(particles[[8, 1]]) @@ -689,7 +822,9 @@ def test17(self): def test18(self): print("SSE validation") - sse_src_path = os.path.join(os.path.dirname(sys.modules[Sse.__module__].__file__), 'src') + sse_src_path = os.path.join( + os.path.dirname(sys.modules[Sse.__module__].__file__), "src" + ) if not os.path.exists(os.path.join(sse_src_path, "evolve.in")): self.skip("Not in a source release") instance = Sse() @@ -701,7 +836,9 @@ def test18(self): instance.stop() testpath = get_path_to_results() - shutil.copy(os.path.join(sse_src_path, "evolve.in"), os.path.join(testpath, "evolve.in")) + shutil.copy( + os.path.join(sse_src_path, "evolve.in"), os.path.join(testpath, "evolve.in") + ) call([os.path.join(sse_src_path, "sse")], cwd=testpath) @@ -709,17 +846,43 @@ def test18(self): lines = sse_output.readlines() sse_final_result = lines[-2].split() - self.assertAlmostEqual(evolved_star.age, float(sse_final_result[0]) | units.Myr, 3) - self.assertAlmostEqual(evolved_star.stellar_type, float(sse_final_result[1]) | units.stellar_type, 3) - self.assertAlmostEqual(evolved_star.initial_mass, float(sse_final_result[2]) | units.MSun, 3) - self.assertAlmostEqual(evolved_star.mass, float(sse_final_result[3]) | units.MSun, 3) - self.assertAlmostEqual(evolved_star.luminosity, 10**float(sse_final_result[4]) | units.LSun, 3) - self.assertAlmostEqual(evolved_star.radius, 10**float(sse_final_result[5]) | units.RSun, 3) - self.assertAlmostRelativeEqual(evolved_star.temperature, 10**float(sse_final_result[6]) | units.K, 2) - self.assertAlmostEqual(evolved_star.core_mass, float(sse_final_result[7]) | units.MSun, 3) - self.assertAlmostEqual(evolved_star.convective_envelope_mass, float(sse_final_result[8]) | units.MSun, 3) - self.assertAlmostEqual(evolved_star.epoch, float(sse_final_result[9]) | units.Myr, 3) - self.assertAlmostEqual(evolved_star.spin, float(sse_final_result[10]) | units.yr**-1, 3) + self.assertAlmostEqual( + evolved_star.age, float(sse_final_result[0]) | units.Myr, 3 + ) + self.assertAlmostEqual( + evolved_star.stellar_type, + float(sse_final_result[1]) | units.stellar_type, + 3, + ) + self.assertAlmostEqual( + evolved_star.initial_mass, float(sse_final_result[2]) | units.MSun, 3 + ) + self.assertAlmostEqual( + evolved_star.mass, float(sse_final_result[3]) | units.MSun, 3 + ) + self.assertAlmostEqual( + evolved_star.luminosity, 10 ** float(sse_final_result[4]) | units.LSun, 3 + ) + self.assertAlmostEqual( + evolved_star.radius, 10 ** float(sse_final_result[5]) | units.RSun, 3 + ) + self.assertAlmostRelativeEqual( + evolved_star.temperature, 10 ** float(sse_final_result[6]) | units.K, 2 + ) + self.assertAlmostEqual( + evolved_star.core_mass, float(sse_final_result[7]) | units.MSun, 3 + ) + self.assertAlmostEqual( + evolved_star.convective_envelope_mass, + float(sse_final_result[8]) | units.MSun, + 3, + ) + self.assertAlmostEqual( + evolved_star.epoch, float(sse_final_result[9]) | units.Myr, 3 + ) + self.assertAlmostEqual( + evolved_star.spin, float(sse_final_result[10]) | units.yr**-1, 3 + ) def test19(self): print("SSE core_mass and CO_core_mass (high mass star)") @@ -760,7 +923,9 @@ def test19(self): def test20(self): print("SSE core_mass and CO_core_mass (low mass stars)") instance = Sse() - stars = instance.particles.add_particles(Particles(mass=[0.6, 1.0] | units.MSun)) + stars = instance.particles.add_particles( + Particles(mass=[0.6, 1.0] | units.MSun) + ) instance.evolve_model(100 | units.Gyr) self.assertEqual(str(stars[0].stellar_type), "Helium White Dwarf") self.assertAlmostEqual(stars[0].mass, 0.405 | units.MSun, 2) @@ -776,17 +941,23 @@ def test21(self): instance = Sse() stars = instance.particles.add_particles(Particles(mass=30 | units.MSun)) mass_loss_wind = stars[0].mass_loss_wind - self.assertAlmostRelativeEquals(mass_loss_wind, 1.703e-07 | units.MSun / units.yr, 3) + self.assertAlmostRelativeEquals( + mass_loss_wind, 1.703e-07 | units.MSun / units.yr, 3 + ) instance.evolve_model(1 | units.Myr) dm = (1 | units.Myr) * mass_loss_wind - self.assertAlmostRelativeEquals(stars[0].mass, (30 | units.MSun) - dm, 3) - self.assertAlmostRelativeEquals(stars[0].mass_loss_wind, 2.053e-07 | units.MSun / units.yr, 3) + self.assertAlmostRelativeEquals(stars[0].mass, (30 | units.MSun) - dm, 3) + self.assertAlmostRelativeEquals( + stars[0].mass_loss_wind, 2.053e-07 | units.MSun / units.yr, 3 + ) instance.stop() def test22(self): instance = Sse() - stars = instance.particles.add_particles(Particles(mass=[1.0, 10.0] | units.MSun)) + stars = instance.particles.add_particles( + Particles(mass=[1.0, 10.0] | units.MSun) + ) gyration_radius = stars.gyration_radius self.assertTrue(numpy.all(0.0 < gyration_radius)) self.assertTrue(numpy.all(gyration_radius < 1.0))