Python interface.EVtwin类代码示例

 def evtwin(self, number_of_stars):
     result = EVtwin()
     result.initialize_code()
     if number_of_stars > result.parameters.maximum_number_of_stars:
         result.parameters.maximum_number_of_stars = number_of_stars
         warnings.warn("You're simulating a large number of stars with EVtwin. This may not be such a good idea...")
     return result

 def xtest19(self):
     print "Testing EVtwin central_temperature and central_density"
     instance = EVtwin()
     stars = instance.particles.add_particles(Particles(mass=[0.1, 1, 10]|units.MSun))
     self.assertIsOfOrder(stars.central_temperature, [4e6, 13e6, 31e6] | units.K)
     self.assertIsOfOrder(stars.central_density, [400, 77, 9] | units.g * units.cm**-3)
     instance.stop()

    def slowtest11(self):
        print "Test for importing new stellar models, also check long-term evolution"
        instance = EVtwin()
        instance.parameters.verbosity = True
        instance.particles.add_particles(Particles(mass = [0.8] | units.MSun))
        instance.evolve_model()

        copy = Particles(1)
        copy.internal_structure = instance.particles[0].get_internal_structure()
        instance.particles.add_particles(copy)

        number_of_zones = instance.particles[0].get_number_of_zones()
        self.assertEqual(len(instance.particles), 2)
        self.assertEqual(instance.particles[1].get_number_of_zones(), number_of_zones)
        self.assertIsOfOrder(instance.particles[1].get_radius_profile()[-1],          1.0 | units.RSun)
        self.assertIsOfOrder(instance.particles[1].get_temperature_profile()[0],  1.0e7 | units.K)
        self.assertIsOfOrder(instance.particles[1].get_pressure_profile()[0],      1.0e17 | units.barye)

        t1, l1 = simulate_evolution_tracks(instance.particles[0], end_time=26.5|units.Gyr)
        t2, l2 = simulate_evolution_tracks(instance.particles[1], end_time=26.5|units.Gyr)
        instance.stop()

        i1 = t1.argmax() # Maximum temperature ~ turnoff main sequence
        i2 = t2.argmax()
        self.assertAlmostRelativeEqual(t1[i1], t2[i2], 2)
        self.assertAlmostRelativeEqual(l1[i1], l2[i2], 2)

 def slowtest7(self):
     print "Test handle_collision"
     position_offset = [100.0, 200.0, 300.0] | units.RSun
     velocity_offset = [10000.0, 20000.0, 30000.0] | units.km / units.s
     colliders = self.new_colliders()
     colliders.position += position_offset
     colliders.velocity += velocity_offset
     
     class GravityCodeStub(object):
         def __init__(self, particles):
             self.particles = particles
     gravity = GravityCodeStub(colliders)
     
     stellar = EVtwin()
     stellar.particles.add_particles(colliders)
     
     collision = StellarEncounterInHydrodynamics(
         #28000, # Bit too slow, even for a slowtest... 
         2800, 
         new_plotting_hydrodynamics_code(Gadget2, 0.2|units.hour, 
             plot_function = pynbody_column_density_plot if HAS_PYNBODY else None,
             plot_function_arguments = dict(width=20|units.RSun, vmin=29, vmax=35) if HAS_PYNBODY else dict(width=20|units.RSun)), 
         hydrodynamics_arguments = dict(redirection="file", redirect_file="hydro_collision_slowtest7_gadget_out.log"),
         hydrodynamics_parameters = dict(time_limit_cpu=1|units.day, gas_epsilon=0.01|units.RSun),
         verbose=True,
         debug=False
     )
     result = collision.handle_collision(colliders[0], colliders[1], gravity_code=gravity, stellar_evolution_code=stellar)
     stellar.stop()
     print result
     self.assertTrue(isinstance(result, Particles))
     self.assertEqual(len(result), 2)
     self.assertTrue((result.mass < [5.0, 2.0] | units.MSun).all())
     self.assertAlmostRelativeEqual(result.center_of_mass(), position_offset, 2)
     self.assertAlmostRelativeEqual(result.center_of_mass_velocity(), velocity_offset, 2)

 def test6(self):
     print "Test handle_collision"
     position_offset = [100.0, 200.0, 300.0] | units.RSun
     velocity_offset = [10000.0, 20000.0, 30000.0] | units.km / units.s
     colliders = self.new_colliders()
     colliders.position *= 1.5 # Grazing collision
     colliders.position += position_offset
     colliders.velocity += velocity_offset
     
     class GravityCodeStub(object):
         def __init__(self, particles):
             self.particles = particles
     gravity = GravityCodeStub(colliders)
     
     stellar = EVtwin()
     stellar.particles.add_particles(colliders)
     
     collision = StellarEncounterInHydrodynamics(280,  # For speed
         Gadget2,
         relax_sph_models=False,  # For speed
         verbose=True)
     collision.extra_steps_when_encounter_is_over = 0 # For speed, no extra steps...
     collision.dynamical_timescales_per_step = 1.3 # ... but then we need to evolve a bit longer in the first step
     result = collision.handle_collision(colliders[0], colliders[1], gravity_code=gravity, stellar_evolution_code=stellar)
     stellar.stop()
     print result
     self.assertTrue(isinstance(result, Particles))
     self.assertEqual(len(result), 2)
     self.assertAlmostEqual(result.mass, [5.0, 2.0] | units.MSun, 1)
     self.assertAlmostRelativeEqual(result.center_of_mass(), position_offset, 2)
     self.assertAlmostRelativeEqual(result.center_of_mass_velocity(), velocity_offset, 2)

 def test5(self):
     print "Test group_bound_particles"
     colliders = self.new_colliders()
     colliders.position = [[0.0, 0.0, 0.0], [1.1, 0.0, 0.0]] | units.RSun
     colliders.velocity = [[0.0, 0.0, 0.0], [10000, 0.0, 0.0]] | units.km / units.s
     stellar = EVtwin()
     stellar.particles.add_particles(colliders)
     
     collision = StellarEncounterInHydrodynamics(7000, None, relax_sph_models=False, verbose=True, 
         star_to_sph_arguments=dict(base_grid_options=dict(type="sobol")))
     collision.dynamical_timescale = zero
     gas_particles = collision.convert_stars(colliders, stellar)
     stellar.stop()
     
     collision.hop = collision.new_hop(gas_particles)
     collision.start_kepler(7 | units.MSun, 10 | units.RSun)
     self.assertTrue(collision.encounter_is_over(gas_particles))
     collision.hop.stop()
     collision.kepler.stop()
     
     groups = collision.groups_after_encounter
     self.assertEqual(len(groups), 2)
     self.assertTrue(4500 < len(groups[0]) < 5000)
     self.assertTrue(1800 < len(groups[1]) < 2000)
     self.assertEqual(len(gas_particles - groups[0] - groups[1]), 346)
     self.assertAlmostEqual(groups[0].center_of_mass()[0], 0 | units.RSun, 1)
     self.assertAlmostEqual(groups[1].center_of_mass()[0], 1.1 | units.RSun, 0)
     self.assertIsOfOrder(groups[1].center_of_mass_velocity()[0], 10000 | units.km / units.s)

    def xtest18(self):
        print "Testing EVtwin calculate_core_mass"
        instance = EVtwin()#redirection="none")
        star = instance.particles.add_particle(Particle(mass=1|units.MSun))
        instance.evolve_model(0.4|units.Gyr) # VERY short, for test speed up
        central_hydrogen_abundance = star.get_chemical_abundance_profiles()[0][0]
        self.assertTrue(central_hydrogen_abundance < 0.68) # some hydrogen is burned
        self.assertTrue(central_hydrogen_abundance > 0.67) # ... but not that much yet
        self.assertEqual(star.calculate_core_mass(core_H_abundance_limit=0.67), 0 | units.MSun)
        self.assertAlmostEqual(star.calculate_core_mass(core_H_abundance_limit=0.71), 1 | units.MSun, 1)

        # For test speed up, we use a weird core_H_abundance_limit to define the "hydrogen exhausted core"
        limit = 0.68
        expected_core_mass = 0.0123182798542 | units.MSun
        self.assertAlmostEqual(star.calculate_core_mass(core_H_abundance_limit=limit), expected_core_mass, 3)

        species_names = star.get_names_of_species()
        self.assertEquals(species_names, ['h1', 'he4', 'c12', 'n14', 'o16', 'ne20', 'mg24', 'si28', 'fe56'])
        h1_core_mass = star.calculate_core_mass(species=["h1"], core_H_abundance_limit=limit)
        he4_core_mass = star.calculate_core_mass(species=["he4"], core_H_abundance_limit=limit)
        c12_core_mass = star.calculate_core_mass(species=["c12"], core_H_abundance_limit=limit)
        n14_core_mass = star.calculate_core_mass(species=["n14"], core_H_abundance_limit=limit)
        o16_core_mass = star.calculate_core_mass(species=["o16"], core_H_abundance_limit=limit)
        ne20_core_mass = star.calculate_core_mass(species=["ne20"], core_H_abundance_limit=limit)
        mg24_core_mass = star.calculate_core_mass(species=["mg24"], core_H_abundance_limit=limit)
        si28_core_mass = star.calculate_core_mass(species=["si28"], core_H_abundance_limit=limit)
        fe56_core_mass = star.calculate_core_mass(species=["fe56"], core_H_abundance_limit=limit)
        metal_core_mass = star.calculate_core_mass(species=["c12", "n14", "o16", "ne20", "mg24", "si28", "fe56"], core_H_abundance_limit=limit)
        instance.stop()
        self.assertAlmostRelativeEqual(h1_core_mass, expected_core_mass*0.68, 2)
        self.assertAlmostRelativeEqual(he4_core_mass, expected_core_mass*0.30, 2)
        self.assertAlmostRelativeEqual(metal_core_mass, expected_core_mass*0.02, 1)
        self.assertAlmostRelativeEqual(expected_core_mass, he4_core_mass + metal_core_mass + h1_core_mass, 7)
        self.assertAlmostRelativeEqual(metal_core_mass, c12_core_mass + n14_core_mass +
            o16_core_mass + ne20_core_mass + mg24_core_mass + si28_core_mass + fe56_core_mass, 7)

def convert_star_to_hydro_model(M, t_end):
    stellar_evolution = EVtwin()
    star = stellar_evolution.particles.add_particle(Particle(mass=M))
    stellar_evolution.evolve_model(t_end)
    Ngas = 10000
    sph_particles = convert_stellar_model_to_SPH(star, Ngas).gas_particles
    stellar_evolution.stop()
    return sph_particles

def set_up_stellar_evolution_code(stars):
    stellar_evolution = EVtwin()
    stellar_evolution.initialize_code()
    # if you run with mesa, you can play with the wind efficiency
    # stellar_evolution.parameters.RGB_wind_scheme = 1
    # stellar_evolution.parameters.reimers_wind_efficiency = 1.0e6 # ridiculous, but instructive
    stellar_evolution.particles.add_particles(stars)
    return stellar_evolution

 def test1(self):
     print "Testing assigned default parameter values..."
     instance = EVtwin()
     instance.initialize_code()
     instance.parameters.set_defaults()
     self.assertEquals(10, instance.parameters.maximum_number_of_stars)
     instance.parameters.maximum_number_of_stars = 12
     self.assertEquals(12, instance.parameters.maximum_number_of_stars)
     instance.stop()

    def test11(self):
        print "Test for importing new stellar models"
        instance = EVtwin()
#~        instance.parameters.import_model_entropy_force = 1.0
#~        instance.parameters.import_model_entropy_accuracy = 1.0e-1
        instance.parameters.verbosity = True
        instance.particles.add_particles(Particles(mass = [0.8] | units.MSun))
        instance.evolve_model()

        instance.new_particle_from_model(instance.particles[0].get_internal_structure())
        # The above line is equivalent with:
        #copy = Particles(1)
        #copy.internal_structure = instance.particles[0].get_internal_structure()
        #instance.particles.add_particles(copy)

        number_of_zones = instance.particles[0].get_number_of_zones()
        self.assertEqual(len(instance.particles), 2)
        self.assertEqual(instance.particles[1].get_number_of_zones(), number_of_zones)
        self.assertIsOfOrder(instance.particles[1].get_radius_profile()[-1],          1.0 | units.RSun)
        self.assertIsOfOrder(instance.particles[1].get_temperature_profile()[0],  1.0e7 | units.K)
        self.assertIsOfOrder(instance.particles[1].get_pressure_profile()[0],      1.0e17 | units.barye)

        instance.evolve_model(keep_synchronous = False)
        self.assertAlmostRelativeEqual(
            instance.particles[0].temperature, instance.particles[1].temperature, 2)
        self.assertAlmostRelativeEqual(
            instance.particles[0].luminosity, instance.particles[1].luminosity, 2)
        instance.stop()

    def xtest13(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 = EVtwin()
        instance.particles.add_particles(stars)

        self.assertAlmostEqual(instance.particles.age, [0.0, 0.0, 0.0] | units.yr)
        self.assertAlmostEqual(instance.particles.time_step, [70465.105509, 6063.68785133, 1876.53255132] | units.yr, 3)

        print "evolve_model without arguments: use shared timestep = 0.99*min(particles.time_step)"
        instance.evolve_model()
        self.assertAlmostEqual(instance.particles.age, 0.99*([1876.53255132,1876.53255132,1876.53255132] | units.yr), 3)
        self.assertAlmostEqual(instance.particles.time_step, [70465.105509,6063.68785133,1876.53255132] | units.yr, 3)
        self.assertAlmostEqual(instance.model_time, 0.99*1876.53255132 | units.yr, 3)

        print "evolve_model with end_time: take timesteps, until end_time is reached exactly"
        instance.evolve_model(15000 | units.yr)
        self.assertAlmostEqual(instance.particles.age, [15000.0, 15000.0, 15000.0] | units.yr, 3)
        self.assertAlmostEqual(instance.particles.time_step, [ 84558.1266108,7276.4254216,2251.83906159] | units.yr, 3)
        self.assertAlmostEqual(instance.model_time, 15000.0 | units.yr, 3)

        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, (15000 | units.yr) + ([ 84558.1266108,7276.4254216,2251.83906159] | units.yr), 3)
        self.assertAlmostRelativeEquals(instance.particles.time_step, [101469.751933,8731.71050591,2702.2068739] | units.yr, 1)
        self.assertAlmostEqual(instance.model_time, 15000.0 | units.yr, 3) # Unchanged!
        instance.stop()

 def test6(self):
     print "Test for obtaining the stellar structure model"
     stars = Particles(2)
     stars.mass = [1.0, 10.0] | units.MSun
     instance = EVtwin()
     instance.initialize_code()
     instance.commit_parameters()
     instance.particles.add_particles(stars)
     instance.commit_particles()
     self.assertEquals(instance.particles.get_number_of_zones(), [199, 199])
     self.assertEquals(len(instance.particles[0].get_radius_profile()), 199)
     self.assertRaises(AmuseException, instance.particles.get_radius_profile,
         expected_message = "Querying radius profiles of more than one particle at a time is not supported.")
     self.assertEquals(len(instance.particles[1].get_density_profile()), 199)
     self.assertIsOfOrder(instance.particles[0].get_radius_profile()[-1],          1.0 | units.RSun)
     self.assertIsOfOrder(instance.particles[0].get_temperature_profile()[0],  1.0e7 | units.K)
     self.assertIsOfOrder(instance.particles[0].get_temperature_profile()[-1],  5.0e3 | units.K)
     radius1 = instance.particles[0].get_radius_profile()
     radius2 = radius1[:-1]
     radius2.prepend(0|units.m)
     delta_radius_cubed = (radius1**3 - radius2**3)
     total_mass = (4./3. * pi * instance.particles[0].get_density_profile() * delta_radius_cubed).sum()
     self.assertAlmostRelativeEqual(total_mass, stars[0].mass, places = 1)
     self.assertAlmostEquals(instance.particles[0].get_mu_profile()[:100], [0.62]*100 | units.amu, places=1)
     instance.stop()

def create_particles():
    star =  Particle()
    star.mass = 3.0 | units.MSun
    
    stellar_evolution = EVtwin()
    se_star = stellar_evolution.particles.add_particle(star)
    
    print("Evolving", star.mass, "star with", stellar_evolution.__class__.__name__, "up to", 100 | units.Myr)
    stellar_evolution.evolve_model(100 | units.Myr)
    
    print("Creating SPH particles from the (1D) stellar evolution model")
    sph_particles = convert_stellar_model_to_SPH(
        se_star, 
        1000
    ).gas_particles
    stellar_evolution.stop()
    return sph_particles

    def slowtest17(self):
        """
        We add multiple particles to evtwin and evolve the stars
        individualy. Evtwin crashes on some combinaties of
        star masses and which star is evolved first.
        """
        exceptions = []
        for i, (masses, indices) in enumerate(
            [
                ([1.21372730283, 1.22207032494, 11.21372730283] | units.MSun, (0, 1)),
                ([1.21372730283, 1.22207032494, 1.21372730283] | units.MSun, (0, 1)),
                ([1.21372730283, 1.22207032494, 1.21372730283] | units.MSun, (1, 0)),
                ([1.21372730283, 1.22207032494] | units.MSun, (0, 1)),
                ([1.21372730283, 11.22207032494, 1.21372730283] | units.MSun, (0, 1)),
                ([1.21372730283, 1.21372730283, 1.21372730283] | units.MSun, (0, 1)),
                ([1.21372730283, 1.22207032494, 1.21372730283] | units.MSun, (0, 1)),
                ([0.101, 1.22207032494, 1.21372730283] | units.MSun, (0, 1)),
            ]
        ):

            stars = Particles(mass=masses)

            instance = EVtwin()
            instance.particles.add_particles(stars)

            try:
                index_in_indices = 0
                instance.particles[indices[index_in_indices]].evolve_for(0.1 | units.Myr)
                index_in_indices = 1
                instance.particles[indices[index_in_indices]].evolve_for(0.1 | units.Myr)
            except AmuseException as ex:
                exceptions.append([i, masses, indices, index_in_indices])

            instance.stop()

        if len(exceptions) > 0:
            failure_message = ""
            for index, masses, indices, index_in_indices in exceptions:
                failure_message += "[{0}]: error in ".format(index)
                failure_message += "first" if index_in_indices == 0 else "second"
                failure_message += " evolve_for,"
                failure_message += " index: {0},".format(indices[index_in_indices])
                failure_message += " masses (MSun): {0}.\n".format(masses.value_in(units.MSun))

            self.fail(failure_message)