Révision 175

NBody/NBody.py (revision 175)
25 25 
    Marsaglia={'CONG':0,'SHR3':1,'MWC':2,'KISS':3}
26 26 
    Computing={'FP32':0,'FP64':1}
27 27 
    Interaction={'Force':0,'Potential':1}
28 
    Artevasion={'None':0,'NegExp':1}
28 
    Artevasion={'None':0,'NegExp':1,'CorRad':2}
29 29 
    return(Marsaglia,Computing,Interaction,Artevasion)
30 30 

  		





  31
  31
  
    
BlobOpenCL= """


  		





  ......




  37
  37
  
    

  		





  38
  38
  
    
#define NONE 0


  		





  39
  39
  
    
#define NEGEXP 1


  		





  
  40
  
    
#define CORRAD 2


  		





  40
  41
  
    

  		





  41
  42
  
    
#if TYPE == TFP32


  		





  42
  43
  
    
#define MYFLOAT4 float4


  		





  ......




  69
  70
  
    

  		





  70
  71
  
    
#define SMALL_NUM (MYFLOAT)1.e-9f


  		





  71
  72
  
    

  		





  
  73
  
    
#define CoreRadius (MYFLOAT)(1.e0f)


  		





  
  74
  
    

  		





  72
  75
  
    
// Create my own Distance implementation: distance buggy on Oland AMD chipset


  		





  73
  76
  
    

  		





  74
  77
  
    
MYFLOAT MyDistance(MYFLOAT4 n,MYFLOAT4 m)


  		





  ......




  83
  86
  
    
    return(sqrt(x2+y2+z2));


  		





  84
  87
  
    
}


  		





  85
  88
  
    

  		





  86
  
  
    
// Interaction to avoid artevasion of bodies : 1-e^-x*x is x*x near 0


  		





  87
  
  
    

  		





  88
  
  
    
//MYFLOAT4 InteractionCore(MYFLOAT4 m,MYFLOAT4 n)


  		





  89
  
  
    
//{


  		





  90
  
  
    
//    private MYFLOAT r=MyDistance((MYFLOAT4)n,(MYFLOAT4)m);


  		





  91
  
  
    
//    private MYFLOAT r2=r*r;


  		





  92
  
  
    
//    private MYFLOAT c=1.e0f/(MYFLOAT)get_global_size(0);


  		





  93
  
  
    

  		





  94
  
  
    
//    private MYFLOAT4 unity=normalize(n,m);


  		





  95
  
  
    

  		





  96
  
  
    
//    return(((MYFLOAT4)n-(MYFLOAT4)m)*(MYFLOAT)(1.e0f-exp(-c*r2))/(MYFLOAT)(r*r2));


  		





  97
  
  
    
//}


  		





  98
  
  
    

  		





  99
  89
  
    
// Potential between 2 m,n bodies


  		





  100
  90
  
    
MYFLOAT PairPotential(MYFLOAT4 m,MYFLOAT4 n)


  		





  101
  91
  
    
#if ARTEVASION == NEGEXP


  		





  ......




  104
  94
  
    
    MYFLOAT r=DISTANCE(n,m);


  		





  105
  95
  
    
    return((-1.e0f+exp(-r))/r);


  		





  106
  96
  
    
}


  		





  
  97
  
    
#elif ARTEVASION == CORRAD


  		





  
  98
  
    
// Add Core Radius to avoid divergence for low distances


  		





  
  99
  
    
{


  		





  
  100
  
    
    MYFLOAT r=DISTANCE(n,m);


  		





  
  101
  
    
    return(-1.e0f/sqrt(r*r+CoreRadius*CoreRadius));


  		





  
  102
  
    
}


  		





  107
  103
  
    
#else


  		





  108
  104
  
    
// Classical potential in 1/r


  		





  109
  105
  
    
{


  		





  ......




  115
  111
  
    
// Interaction based of Force (1/r**2) or Potential (-grad(1/r))


  		





  116
  112
  
    
MYFLOAT4 Interaction(MYFLOAT4 m,MYFLOAT4 n)


  		





  117
  113
  
    
#if INTERACTION == TFORCE


  		





  
  114
  
    
#if ARTEVASION == NEGEXP


  		





  
  115
  
    
// Force gradient of potential, set as (1-exp(-r))/r 


  		





  
  116
  
    
{


  		





  
  117
  
    
    private MYFLOAT r=MyDistance(n,m);


  		





  
  118
  
    
    private MYFLOAT num=1.e0f+exp(-r)*(r-1.e0f);


  		





  
  119
  
    
    return((n-m)*num/(MYFLOAT)(r*r*r));


  		





  
  120
  
    
}


  		





  
  121
  
    
#elif ARTEVASION == CORRAD


  		





  
  122
  
    
// Force gradient of potential, (Core Radius) set as 1/sqrt(r**2+CoreRadius**2) 


  		





  
  123
  
    
{


  		





  
  124
  
    
    private MYFLOAT r=MyDistance(n,m);


  		





  
  125
  
    
    private MYFLOAT den=sqrt(r*r+CoreRadius*CoreRadius);


  		





  
  126
  
    
    return((n-m)/(MYFLOAT)(den*den*den));


  		





  
  127
  
    
}


  		





  
  128
  
    
#else


  		





  118
  129
  
    
// Simplest implementation of force (equals to acceleration)


  		





  119
  130
  
    
// seems to bo bad (numerous artevasions)


  		





  120
  131
  
    
// MYFLOAT4 InteractionForce(MYFLOAT4 m,MYFLOAT4 n)


  		





  ......




  122
  133
  
    
    private MYFLOAT r=MyDistance(n,m);


  		





  123
  134
  
    
    return((n-m)/(MYFLOAT)(r*r*r));


  		





  124
  135
  
    
}


  		





  
  136
  
    
#endif


  		





  125
  137
  
    
#else


  		





  126
  138
  
    
// Force definited as gradient of potential


  		





  127
  139
  
    
// Estimate potential and proximate potential to estimate force


  		





  ......




  523
  535
  
    
    Switch = as_8_bit( 's' )


  		





  524
  536
  
    

  		





  525
  537
  
    
    Zoom=1


  		





  526
  
  
    
    if k == GLUT_KEY_UP:


  		





  527
  
  
    
        ViewRX += 1.0


  		





  528
  
  
    
    elif k == GLUT_KEY_DOWN:


  		





  529
  
  
    
        ViewRX -= 1.0


  		





  530
  
  
    
    elif k == GLUT_KEY_LEFT:


  		





  531
  
  
    
        ViewRY += 1.0


  		





  532
  
  
    
    elif k == GLUT_KEY_RIGHT:


  		





  533
  
  
    
        ViewRY -= 1.0


  		





  534
  
  
    
    elif k == LC_Z:


  		





  
  538
  
    
    if k == LC_Z:


  		





  535
  539
  
    
        ViewRZ += 1.0


  		





  536
  540
  
    
    elif k == UC_Z:


  		





  537
  541
  
    
        ViewRZ -= 1.0


  		





  ......




  624
  628
  
    
    # Type of Interaction


  		





  625
  629
  
    
    InterType='Force'


  		





  626
  630
  
    
    


  		





  627
  
  
    
    HowToUse='%s -h [Help] -r [InitialRandom] -g [OpenGL] -e [VirielStress] -o [Verbose] -p [Potential] -x [NegativeExponential4ArteEvasions] -d <DeviceId> -n <NumberOfParticules> -i <Iterations> -z <SizeOfBoxOrBall> -v <Velocity> -s <Step> -b <Ball|Box> -m <ImplicitEuler|RungeKutta|ExplicitEuler|Heun> -t <FP32|FP64>'


  		





  
  631
  
    
    HowToUse='%s -h [Help] -r [InitialRandom] -g [OpenGL] -e [VirielStress] -o [Verbose] -p [Potential] -x <None|NegExp|CorRad> -d <DeviceId> -n <NumberOfParticules> -i <Iterations> -z <SizeOfBoxOrBall> -v <Velocity> -s <Step> -b <Ball|Box> -m <ImplicitEuler|RungeKutta|ExplicitEuler|Heun> -t <FP32|FP64>'


  		





  628
  632
  
    

  		





  629
  633
  
    
    try:


  		





  630
  
  
    
        opts, args = getopt.getopt(sys.argv[1:],"rpxgehod:n:i:z:v:s:m:t:b:",["random","potential","negexp","opengl","viriel","verbose","device=","number=","iterations=","size=","velocity=","step=","method=","valuetype=","shape="])


  		





  
  634
  
    
        opts, args = getopt.getopt(sys.argv[1:],"rpgehod:n:i:z:v:s:m:t:b:x:",["random","potential","coarev","opengl","viriel","verbose","device=","number=","iterations=","size=","velocity=","step=","method=","valuetype=","shape="])


  		





  631
  635
  
    
    except getopt.GetoptError:


  		





  632
  636
  
    
        print(HowToUse % sys.argv[0])


  		





  633
  637
  
    
        sys.exit(2)


  		





  ......




  663
  667
  
    
            Method=arg


  		





  664
  668
  
    
        elif opt in ("-b", "--shape"):


  		





  665
  669
  
    
            Shape=arg


  		





  
  670
  
    
            if Shape!='Ball' or Shape!='Box':


  		





  
  671
  
    
                print('Wrong argument: set to Ball')


  		





  666
  672
  
    
        elif opt in ("-n", "--number"):


  		





  667
  673
  
    
            Number=int(arg)


  		





  668
  674
  
    
        elif opt in ("-i", "--iterations"):


  		





  ......




  684
  690
  
    
            OpenGL=True


  		





  685
  691
  
    
        elif opt in ("-p", "--potential"):


  		





  686
  692
  
    
            InterType='Potential'


  		





  687
  
  
    
        elif opt in ("-x", "--negexp"):


  		





  688
  
  
    
            CoArEv='NegExp'


  		





  
  693
  
    
        elif opt in ("-x", "--coarev"):


  		





  
  694
  
    
            CoArEv=arg


  		





  689
  695
  
    
        elif opt in ("-o", "--verbose"):


  		





  690
  696
  
    
            Verbose=True


  		





  691
  697
  
    

  		





  692
  
  
    
    SizeOfShape=MyFloat(SizeOfShape*Number)


  		





  
  698
  
    
    SizeOfShape=np.sqrt(MyFloat(SizeOfShape*Number))


  		





  693
  699
  
    
    Velocity=MyFloat(Velocity)


  		





  694
  700
  
    
    Step=MyFloat(Step)


  		





  695
  701
  
    
                


  		












Formats disponibles :
  Unified diff