Python interface.Huayno类代码示例

    def test3(self):
        convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 149.5e6 | units.km)

        instance = Huayno(convert_nbody)
        instance.initialize_code()
        instance.parameters.epsilon_squared = 0.00001 | units.AU**2
        
        stars = datamodel.Stars(2)
        star1 = stars[0]
        star2 = stars[1]

        star1.mass = units.MSun(1.0)
        star1.position = units.AU(numpy.array((-1.0,0.0,0.0)))
        star1.velocity = units.AUd(numpy.array((0.0,0.0,0.0)))
        star1.radius = units.RSun(1.0)

        star2.mass = units.MSun(1.0)
        star2.position = units.AU(numpy.array((1.0,0.0,0.0)))
        star2.velocity = units.AUd(numpy.array((0.0,0.0,0.0)))
        star2.radius = units.RSun(100.0)
        
        instance.particles.add_particles(stars)
    
        for x in range(1,2000,10):
            instance.evolve_model(x | units.day)
            instance.particles.copy_values_of_all_attributes_to(stars)
            stars.savepoint()

 def test16(self):
     instance = Huayno()
     instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2
     
     particles = datamodel.Particles(2)
     particles.mass = [1.] | nbody_system.mass
     particles.radius =  [0.0] | nbody_system.length
     particles.position = [[0.0,0.0,0.0],[1.0,0.0,0.0]] | nbody_system.length
     particles.velocity = [[0.0, 0.0, 0.0]] | nbody_system.speed
     instance.particles.add_particles(particles)
     
     zero = 0.0 | nbody_system.length
     
     for x in (0.25, 0.5, 0.75):
         x0 = x | nbody_system.length
         potential0 = instance.get_potential_at_point(zero, x0, zero, zero)
         fx0, fy0, fz0 = instance.get_gravity_at_point(zero, x0, zero, zero)
         
         self.assertAlmostEqual(fy0, 0.0 | nbody_system.acceleration,14)
         self.assertAlmostEqual(fz0, 0.0 | nbody_system.acceleration,14)
         
         fx = (-1.0 / (x0**2)+1.0 / (((1.0|nbody_system.length)-x0)**2)) * (1.0 | nbody_system.length ** 3 / nbody_system.time ** 2)
         self.assertAlmostEqual(fx, fx0,14)
         self.assertAlmostEqual(potential0, -nbody_system.G*(1.|nbody_system.mass)*(1./x0+1./((1.|nbody_system.length)-x0)),14)
     instance.stop()

    def test16(self):
        solsys = new_solar_system()

        solsys.x-=1.| units.AU

        p=datamodel.Particles(3)
        p.mass=[1,2,3] | units.MSun
        p.x=[1,10,100] | units.AU
        p.y=[0,0,-10] | units.AU
        p.z=[0,0,10] | units.AU

        from amuse.community.huayno.interface import Huayno
        from amuse.units import nbody_system        
        conv=nbody_system.nbody_to_si(1. | units.AU, 1.| units.MSun)
        h = Huayno(conv)
        h.particles.add_particles(solsys)
        
        ax1,ay1,az1=h.get_gravity_at_point(p.x*0.,p.x,p.y,p.z)
        
        mercury = Mercury()        
        mercury.particles.add_particles(solsys)

        ax2,ay2,az2=mercury.get_gravity_at_point(p.x*0.,p.x,p.y,p.z)


        self.assertAlmostRelativeEqual(ax1,ax2,12)
        self.assertAlmostRelativeEqual(ay1,ay2,12)
        self.assertAlmostRelativeEqual(az1,az2,12)

def initialize_huayno(bodies, converter, huayno_eta):
  stars_gravity = Huayno(converter,channel_type="sockets")
  stars_gravity.particles.add_particles(bodies) # either bodies or just stars
  stars_gravity.commit_particles()
  #stars_gravity.set_timestep_parameter(huayno_eta)
  stars_gravity.parameters.timestep_parameter = huayno_eta
  stars_gravity.set_inttype_parameter(12) # CC_KEPLER

  return stars_gravity

    def test5(self):
        convert_nbody = nbody_system.nbody_to_si(5.0 | units.kg, 10.0 | units.m)

        instance = Huayno(convert_nbody)
        instance.initialize_code()
        
        particles = datamodel.Particles(2)
        self.assertEquals(len(instance.particles), 0)
        
        particles.mass = [15.0, 30.0] | units.kg
        particles.radius =  [10.0, 20.0] | units.m
        particles.position = [[10.0, 20.0, 30.0], [20.0, 40.0, 60.0]] | units.m
        particles.velocity = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] | units.m / units.s

        
        instance.particles.add_particles(particles)
        self.assertEquals(len(instance.particles), 2)
        
        instance.set_state(1, 16|units.kg, 20.0|units.m, 40.0|units.m, 60.0|units.m, 
                                 1.0|units.ms, 1.0|units.ms, 1.0|units.ms)
        
        curr_state =  instance.get_state(1)
        for expected, actural in zip((16|units.kg, 20.0|units.m, 40.0|units.m, 60.0|units.m, 
                                 1.0|units.ms, 1.0|units.ms, 1.0|units.ms, 0 | units.m), curr_state):
            self.assertAlmostRelativeEquals(actural,expected)
        
        instance.set_state(1, 16|units.kg, 20.0|units.m, 40.0|units.m, 60.0|units.m, 
                                 1.0|units.ms, 1.0|units.ms, 1.0|units.ms , 20.0|units.m)
        
        curr_state =  instance.get_state(1)
        for expected, actural in zip((16|units.kg, 20.0|units.m, 40.0|units.m, 60.0|units.m, 
                                 1.0|units.ms, 1.0|units.ms, 1.0|units.ms, 20 | units.m), curr_state):
            self.assertAlmostRelativeEquals(actural,expected)

def binary_with_planet_run(m1=1.|units.MSun, m2=1.| units.MSun, m_planet=1|units.MJupiter,
                           r1=None, r2=None, r_planet=None, ecc_binary=0, P_binary=20 | units.day,
                           ecc_planet=0., P_planet=1.| units.yr, pangle_planet=0., a_planet=None,
                           tend=100. | units.yr,hostname=''):

    dEcrit=1.e-6

    three=binary_with_planet(m1=m1,m2=m2,m_planet=m_planet,r1=r1,r2=r2,r_planet=r_planet,
             ecc_binary=ecc_binary,P_binary=P_binary,
             ecc_planet=ecc_planet,a_planet=a_planet,pangle_planet=pangle_planet)
    
    convert=nbody_system.nbody_to_si(1|units.MSun,1|units.AU)
    code=Huayno(convert,hostname=hostname)
    
    code.parameters.inttype_parameter=code.inttypes.SHARED10
    code.parameters.timestep_parameter=.2
    code.parameters.timestep=100. | units.day
        
    dt=10000. | units.day
    
    code.particles.add_particles(three)

    E0=code.potential_energy+code.kinetic_energy
    a0,eps0,pangle0=elements( three.total_mass(),
          code.particles.x[2],
          code.particles.y[2],
          code.particles.z[2],
          code.particles.vx[2],
          code.particles.vy[2],
          code.particles.vz[2] )

    t=0. | units.day
    result="stable"
    while(t < tend-dt/2):
        t=t+dt
        code.evolve_model(t)
        E=code.potential_energy+code.kinetic_energy    
        dE=abs(((E-E0)/E0))
        a,eps,pangle=elements( three.total_mass(),
          code.particles.x[2],
          code.particles.y[2],
          code.particles.z[2],
          code.particles.vx[2],
          code.particles.vy[2],
          code.particles.vz[2] )
        if dE > dEcrit or a<0.5*a0 or a>2.*a0:
          result="unstable"
          if dE > dEcrit:
            result="failed"
          break
    code.stop()
    return result,t,dE,a.in_(units.AU),eps,pangle

    def test14(self):
        import hashlib

        numpy.random.seed(123456)
        particles = plummer.new_plummer_model(64)
        sha=hashlib.sha1()

        class inttypes(object):
            SHARED2=1
            EXTRAPOLATE=5
            PASS_KDK=2
            PASS_DKD=7
            HOLD_KDK=3
            HOLD_DKD=8
            PPASS_DKD=9
            BRIDGE_KDK=4
            BRIDGE_DKD=10
            CC=11
            CC_KEPLER=12
            OK=13
            KEPLER=14
            SHARED4=15
            SHARED6=18
            SHARED8=19
            SHARED10=20
            SHAREDBS=21
           
            @classmethod
            def _list(cls):
                  return set([x for x in cls.__dict__.keys() if not x.startswith('_')])

        for itype in sorted(inttypes._list()):
            if itype in ("KEPLER"): continue
            instance = Huayno()
            instance.parameters.inttype_parameter=getattr(Huayno.inttypes,itype)
            instance.particles.add_particles(particles)
            instance.evolve_model(0.125 | nbody_system.time)
            part_out= instance.particles.copy()
            position = part_out.position.number
            if hasattr(position,'tobytes'):
                as_bytes = position.tobytes()
            else:
                as_bytes = numpy.copy(position.data, order='C')
            sha.update(as_bytes)
            
            instance.stop()
        
        # this result is probably dependent on system architecture hence no good for assert
        print 
        print sha.hexdigest()
        print "8e71f9441578a43af4af927943577ad2c4130e4c"

 def test13(self):
     particles = plummer.new_plummer_model(31)
    
     instance = Huayno(number_of_workers=1)#, debugger="xterm")
     instance.initialize_code()
     instance.parameters.epsilon_squared = 0.01 | nbody_system.length ** 2
     instance.particles.add_particles(particles)
     
     instance.evolve_model(0.1 | nbody_system.time)
     instance.synchronize_model()
     expected_positions = instance.particles.position
     instance.stop()

 def test29(self):
     instance = Huayno()
     instance.parameters.epsilon_squared = 0.0 | nbody_system.length**2
     
     particles = datamodel.Particles(2)
     particles.position = ([0,0,0], [1,0,0] )| nbody_system.length
     particles.velocity = ([-1,0,0], [2,0,0] )| nbody_system.speed
     particles.radius = 0| nbody_system.length
     particles.mass = 0.1| nbody_system.mass
             
     instance.particles.add_particles(particles)
     instance.stopping_conditions.out_of_box_detection.enable()
     instance.parameters.stopping_conditions_out_of_box_size = 2 | nbody_system.length
     instance.parameters.stopping_conditions_out_of_box_use_center_of_mass = False
     instance.evolve_model(1 | nbody_system.time)
     self.assertTrue(instance.stopping_conditions.out_of_box_detection.is_set())
     self.assertEquals(len(instance.stopping_conditions.out_of_box_detection.particles(0)), 1)
     self.assertEquals(instance.stopping_conditions.out_of_box_detection.particles(0)[0].key, particles[1].key)
     instance.stop()

 def _compare_integrator_with_collision_integrator(self, inttype_parameter1, inttype_parameter2):
     numpy.random.seed(12345)
     particles = plummer.new_plummer_model(101)
     instance = Huayno()
     instance.parameters.inttype_parameter = inttype_parameter1
     instance.particles.add_particles(particles)
     instance.evolve_model(0.2 | nbody_system.time)
     expected_position = instance.particles.position
     expected_velocity = instance.particles.velocity
     instance.reset()
     instance.parameters.inttype_parameter = inttype_parameter2
     instance.particles.add_particles(particles)
     instance.evolve_model(0.2 | nbody_system.time)
     self.assertAlmostRelativeEquals(expected_position, instance.particles.position, 8)
     self.assertAlmostRelativeEquals(expected_velocity, instance.particles.velocity, 8)
     instance.stop()

    def test24(self):
        print "test massless particles/ kepler integrator"
        N=20        
        tend=2.| units.yr
        numpy.random.seed(12345)
        conv=nbody_system.nbody_to_si(4.| units.MSun, 5.|units.AU)
        orbiters=plummer.new_plummer_model(N,conv)
        sun=datamodel.Particle(mass=1.|units.MSun)
        sun.position=[12.3,1.,-.2]|units.AU
        sun.velocity=[10,50.,-20.]|units.kms
        orbiters.mass*=0.
        
        a0,eps0=elements(sun.mass,orbiters.x,orbiters.y,orbiters.z,
                      orbiters.vx,orbiters.vy,orbiters.vz)

        orbiters.position+=sun.position
        orbiters.velocity+=sun.velocity

        pos=dict()
        for inttype in [20,14]:
          code=Huayno(conv)
          code.parameters.inttype_parameter=inttype
          code.parameters.timestep_parameter=0.1
          code.particles.add_particle(sun)
          orbiters2=code.particles.add_particles(orbiters).copy()
          orbiters2.position-=sun.position
          orbiters2.velocity-=sun.velocity
          code.evolve_model(tend)
          a,eps=elements(sun.mass,orbiters2.x,orbiters2.y,orbiters2.z,
                    orbiters2.vx,orbiters2.vy,orbiters2.vz)

          da=abs((a-a0)/a0)
          deps=abs(eps-eps0)/eps0

          dev=numpy.where(da > 1.e-12)[0]
          self.assertEqual( len(dev),0)
          dev=numpy.where(deps > 1.e-12)[0]
          self.assertEqual( len(dev),0)
          pos[inttype]=[orbiters2.x.value_in(units.AU),orbiters2.y.value_in(units.AU),orbiters2.z.value_in(units.AU)]
        self.assertAlmostEqual(pos[20][0],pos[14][0],12)
        self.assertAlmostEqual(pos[20][1],pos[14][1],12)
        self.assertAlmostEqual(pos[20][2],pos[14][2],12)

 def test6(self):
     print "Test6: Testing Huayno parameters"
     convert_nbody = nbody_system.nbody_to_si(1.0 | units.yr, 1.0 | units.AU)
     instance = Huayno(convert_nbody)
     
     (value, error) = instance.legacy_interface.get_eps2()
     self.assertEquals(0, error)
     self.assertEquals(0.0, value)
     self.assertAlmostEquals(0.0 | units.AU**2, instance.parameters.epsilon_squared, in_units=units.AU**2)
     for x in [0.01, 0.1, 0.2]:
         instance.parameters.epsilon_squared = x | units.AU**2
         self.assertAlmostEquals(x | units.AU**2, instance.parameters.epsilon_squared, in_units=units.AU**2)
             
     (value, error) = instance.legacy_interface.get_time()
     self.assertEquals(0, error)
     self.assertEquals(0.0, value)
     self.assertAlmostEquals(0.0 | units.yr, instance.parameters.begin_time, in_units=units.yr)
     for x in [1.0, 10.0, 100.0]:
         instance.parameters.begin_time = x | units.yr
         self.assertAlmostEquals(x | units.yr, instance.parameters.begin_time, in_units=units.yr)
     instance.stop()

    def test27(self):
        particles = plummer.new_plummer_model(31)

        tend=0.25| nbody_system.time

        instance = Huayno()
        instance.particles.add_particles(particles)        
        instance.evolve_model(tend)
        expected_positions = instance.particles.position
        self.assertEqual(instance.model_time, tend)
        instance.stop()
        
        particles2=particles.copy()
        particles2.velocity*=-1

        instance = Huayno()
        instance.particles.add_particles(particles2)        
        instance.evolve_model(-tend)
        positions = instance.particles.position
        self.assertEqual(instance.model_time, -tend)
        instance.stop()

        self.assertAlmostEqual(positions,expected_positions)

    def test4(self):
        convert_nbody = nbody_system.nbody_to_si(5.0 | units.kg, 10.0 | units.m)

        instance = Huayno(convert_nbody)
        instance.initialize_code()
        
        particles = datamodel.Particles(2)
        self.assertEquals(len(instance.particles), 0)
        
        particles.mass = [15.0, 30.0] | units.kg
        particles.radius =  [10.0, 20.0] | units.m
        particles.position = [[10.0, 20.0, 30.0], [20.0, 40.0, 60.0]] | units.m
        particles.velocity = [[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]] | units.m / units.s

        
        instance.particles.add_particles(particles)
        self.assertEquals(len(instance.particles), 2)
        
        instance.particles.mass =  [17.0, 33.0] | units.kg
        
        
        self.assertEquals(instance.get_mass(0), 17.0| units.kg) 
        self.assertEquals(instance.get_mass(1), 33.0| units.kg)  

 def test22(self):
     print "Testing zero-mass test particles in Huayno, can be used for removing particles when inside recursive evolve loop"
     sun_and_earth = self.new_system_of_sun_and_earth()
     period = (4.0 * math.pi**2 * (1.0 | units.AU)**3 / (constants.G * sun_and_earth.total_mass())).sqrt()
     convert_nbody = nbody_system.nbody_to_si(1.0 | units.MSun, 1.0 | units.AU)
     huayno = Huayno(convert_nbody)
     huayno.parameters.epsilon_squared = 0.0 | units.AU**2
     huayno.parameters.inttype_parameter = huayno.inttypes.SHARED8
     
     test_particle = datamodel.Particle(mass=0|units.MSun, position=[4,0,0]|units.AU, velocity=[0,0,0]|units.kms)
     test_particle.vy = (constants.G * sun_and_earth.total_mass() / (4.0 | units.AU)).sqrt()
     sun_and_earth.add_particle(test_particle)
     huayno.particles.add_particles(sun_and_earth)
     huayno.evolve_model(period)
     self.assertAlmostRelativeEqual(huayno.particles[:2].x, sun_and_earth[:2].x, 13)
     huayno.evolve_model(1.25 * period)
     self.assertAlmostRelativeEqual(huayno.particles[1].y, sun_and_earth[1].x, 13)
     huayno.evolve_model(8.0 * period)
     self.assertAlmostRelativeEqual(huayno.particles.x, sun_and_earth.x, 8)
     huayno.stop()