Révision 218 TrouNoir/TrouNoir.py
| TrouNoir.py (revision 218) | ||
|---|---|---|
| 227 | 227 |
raie=Line; |
| 228 | 228 |
|
| 229 | 229 |
private float d,bmx,db,b,h; |
| 230 |
private float rp[TRACKPOINTS];
|
|
| 230 |
private float rp0,rpp,rps;
|
|
| 231 | 231 |
private float phi,thi,phd,php,nr,r; |
| 232 | 232 |
private int nh; |
| 233 | 233 |
private float zp,fp; |
| ... | ... | |
| 250 | 250 |
// impact parameter |
| 251 | 251 |
b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx; |
| 252 | 252 |
// step of impact parameter; |
| 253 |
// db=bmx/(float)(sizex/2); |
|
| 254 | 253 |
db=bmx/(float)(sizex); |
| 255 | 254 |
|
| 256 | 255 |
up=0.; |
| ... | ... | |
| 261 | 260 |
|
| 262 | 261 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
| 263 | 262 |
|
| 264 |
rp[nh]=fabs(b/us); |
|
| 263 |
rps=fabs(b/us); |
|
| 264 |
rp0=rps; |
|
| 265 | 265 |
|
| 266 | 266 |
int ExitOnImpact=0; |
| 267 | 267 |
|
| ... | ... | |
| 271 | 271 |
pp=ps; |
| 272 | 272 |
up=us; |
| 273 | 273 |
vp=vs; |
| 274 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 275 |
rp[nh]=fabs(b/us); |
|
| 276 |
ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0; |
|
| 274 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 275 |
rpp=rps; |
|
| 276 |
rps=fabs(b/us); |
|
| 277 |
ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rps>ri)&&(rps<re)?1:0; |
|
| 277 | 278 |
|
| 278 |
} while ((rp[nh]>=rs)&&(rp[nh]<=rp[0])&&(ExitOnImpact==0));
|
|
| 279 |
} while ((rps>=rs)&&(rps<=rp0)&&(ExitOnImpact==0));
|
|
| 279 | 280 |
|
| 280 | 281 |
if (ExitOnImpact==1) {
|
| 281 |
impact(phi,rp[nh-1],b,tho,m,&zp,&fp,q,db,h,raie);
|
|
| 282 |
impact(phi,rpp,b,tho,m,&zp,&fp,q,db,h,raie);
|
|
| 282 | 283 |
} |
| 283 | 284 |
else |
| 284 | 285 |
{
|
| ... | ... | |
| 712 | 713 |
float ExternalRadius,float Angle, |
| 713 | 714 |
int Line) |
| 714 | 715 |
{
|
| 715 |
uint xi=(uint)blockIdx.x;
|
|
| 716 |
uint yi=(uint)blockIdx.y;
|
|
| 716 |
uint xi=(uint)(blockIdx.x*blockDim.x+threadIdx.x);
|
|
| 717 |
uint yi=(uint)(blockIdx.y*blockDim.y+threadIdx.y);
|
|
| 717 | 718 |
uint sizex=(uint)gridDim.x*blockDim.x; |
| 718 | 719 |
uint sizey=(uint)gridDim.y*blockDim.y; |
| 719 | 720 |
|
| ... | ... | |
| 731 | 732 |
raie=Line; |
| 732 | 733 |
|
| 733 | 734 |
float d,bmx,db,b,h; |
| 734 |
float rp[TRACKPOINTS];
|
|
| 735 |
float rp0,rpp,rps;
|
|
| 735 | 736 |
float phi,thi,phd,php,nr,r; |
| 736 | 737 |
int nh; |
| 737 | 738 |
float zp,fp; |
| ... | ... | |
| 765 | 766 |
|
| 766 | 767 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
| 767 | 768 |
|
| 768 |
rp[nh]=fabs(b/us); |
|
| 769 |
rps=fabs(b/us); |
|
| 770 |
rp0=rps; |
|
| 769 | 771 |
|
| 770 | 772 |
int ExitOnImpact=0; |
| 771 | 773 |
|
| ... | ... | |
| 775 | 777 |
pp=ps; |
| 776 | 778 |
up=us; |
| 777 | 779 |
vp=vs; |
| 778 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 779 |
rp[nh]=fabs(b/us); |
|
| 780 |
ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0; |
|
| 780 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 781 |
rpp=rps; |
|
| 782 |
rps=fabs(b/us); |
|
| 783 |
ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rps>ri)&&(rps<re)?1:0; |
|
| 781 | 784 |
|
| 782 |
} while ((rp[nh]>=rs)&&(rp[nh]<=rp[0])&&(ExitOnImpact==0));
|
|
| 785 |
} while ((rps>=rs)&&(rps<=rp0)&&(ExitOnImpact==0));
|
|
| 783 | 786 |
|
| 784 | 787 |
if (ExitOnImpact==1) {
|
| 785 |
impact(phi,rp[nh-1],b,tho,m,&zp,&fp,q,db,h,raie);
|
|
| 788 |
impact(phi,rpp,b,tho,m,&zp,&fp,q,db,h,raie);
|
|
| 786 | 789 |
} |
| 787 | 790 |
else |
| 788 | 791 |
{
|
| ... | ... | |
| 790 | 793 |
fp=0.; |
| 791 | 794 |
} |
| 792 | 795 |
|
| 796 |
__syncthreads(); |
|
| 797 |
|
|
| 793 | 798 |
zImage[yi+sizex*xi]=(float)zp; |
| 794 | 799 |
fImage[yi+sizex*xi]=(float)fp; |
| 795 | 800 |
} |
| ... | ... | |
| 1165 | 1170 |
numpy.int32(Line)) |
| 1166 | 1171 |
CLLaunch.wait() |
| 1167 | 1172 |
|
| 1168 |
elapsed = time.time()-start_time |
|
| 1169 |
print("\nElapsed Time : %f" % elapsed)
|
|
| 1173 |
compute = time.time()-start_time |
|
| 1170 | 1174 |
|
| 1171 | 1175 |
cl.enqueue_copy(queue,zImage,zImageCL).wait() |
| 1172 | 1176 |
cl.enqueue_copy(queue,fImage,fImageCL).wait() |
| 1173 | 1177 |
cl.enqueue_copy(queue,Trajectories,TrajectoriesCL).wait() |
| 1174 | 1178 |
cl.enqueue_copy(queue,IdLast,IdLastCL).wait() |
| 1179 |
elapsed = time.time()-start_time |
|
| 1180 |
print("\nCompute Time : %f" % compute)
|
|
| 1181 |
print("Elapsed Time : %f\n" % elapsed)
|
|
| 1175 | 1182 |
|
| 1176 | 1183 |
zMaxPosition=numpy.where(zImage[:,:]==zImage.max()) |
| 1177 | 1184 |
fMaxPosition=numpy.where(fImage[:,:]==fImage.max()) |
| ... | ... | |
| 1261 | 1268 |
numpy.float32(ExternalRadius), |
| 1262 | 1269 |
numpy.float32(Angle), |
| 1263 | 1270 |
numpy.int32(Line), |
| 1264 |
grid=(zImage.shape[0],zImage.shape[1]),block=(1,1,1))
|
|
| 1271 |
grid=(zImage.shape[0]/32,zImage.shape[1]/32),block=(32,32,1))
|
|
| 1265 | 1272 |
elif Method=='TrajectoCircle': |
| 1266 | 1273 |
TrajectoryCU(TrajectoriesCU,IdLastCU, |
| 1267 | 1274 |
numpy.uint32(Trajectories.shape[1]), |
| ... | ... | |
| 1301 | 1308 |
grid=(zImage.shape[0],zImage.shape[1],1),block=(1,1,1)) |
| 1302 | 1309 |
|
| 1303 | 1310 |
|
| 1304 |
elapsed = time.time()-start_time |
|
| 1305 |
print("\nElapsed Time : %f" % elapsed)
|
|
| 1311 |
compute = time.time()-start_time |
|
| 1306 | 1312 |
|
| 1307 | 1313 |
cuda.memcpy_dtoh(zImage,zImageCU) |
| 1308 | 1314 |
cuda.memcpy_dtoh(fImage,fImageCU) |
| 1315 |
elapsed = time.time()-start_time |
|
| 1316 |
print("\nCompute Time : %f" % compute)
|
|
| 1317 |
print("Elapsed Time : %f\n" % elapsed)
|
|
| 1309 | 1318 |
|
| 1310 | 1319 |
zMaxPosition=numpy.where(zImage[:,:]==zImage.max()) |
| 1311 | 1320 |
fMaxPosition=numpy.where(fImage[:,:]==fImage.max()) |
| 1312 | 1321 |
print("Z max @(%i,%i) : %f" % (zMaxPosition[0][0],zMaxPosition[1][0],zImage.max()))
|
| 1313 | 1322 |
print("Flux max @(%i,%i) : %f\n" % (fMaxPosition[0][0],fMaxPosition[1][0],fImage.max()))
|
| 1314 | 1323 |
|
| 1324 |
Context.pop() |
|
| 1325 |
|
|
| 1315 | 1326 |
Context.detach() |
| 1316 | 1327 |
|
| 1317 | 1328 |
return(elapsed) |
Formats disponibles : Unified diff