Révision 225 TrouNoir/TrouNoir.py
| TrouNoir.py (revision 225) | ||
|---|---|---|
| 532 | 532 |
|
| 533 | 533 |
} |
| 534 | 534 |
|
| 535 |
__kernel void Original(__global float *zImage,__global float *fImage,
|
|
| 536 |
float Mass,float InternalRadius, |
|
| 537 |
float ExternalRadius,float Angle, |
|
| 538 |
int Line) |
|
| 535 |
__kernel void EachCircle(__global float *zImage,__global float *fImage,
|
|
| 536 |
float Mass,float InternalRadius,
|
|
| 537 |
float ExternalRadius,float Angle,
|
|
| 538 |
int Line)
|
|
| 539 | 539 |
{
|
| 540 | 540 |
// Integer Impact Parameter ID |
| 541 | 541 |
int bi=get_global_id(0); |
| 542 | 542 |
// Integer Impact Parameter Size (half of image) |
| 543 | 543 |
int bmaxi=get_global_size(0); |
| 544 | 544 |
|
| 545 |
float Trajectories[2048];
|
|
| 545 |
float Trajectory[TRACKPOINTS];
|
|
| 546 | 546 |
|
| 547 |
// Perform trajectory for each pixel |
|
| 548 |
|
|
| 549 | 547 |
float m,rs,ri,re,tho; |
| 550 | 548 |
int raie,q; |
| 551 | 549 |
|
| ... | ... | |
| 583 | 581 |
// b versus us |
| 584 | 582 |
float bvus=fabs(b/us); |
| 585 | 583 |
float bvus0=bvus; |
| 586 |
Trajectories[nh]=bvus;
|
|
| 584 |
Trajectory[nh]=bvus;
|
|
| 587 | 585 |
|
| 588 | 586 |
do |
| 589 | 587 |
{
|
| ... | ... | |
| 593 | 591 |
vp=vs; |
| 594 | 592 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
| 595 | 593 |
bvus=fabs(b/us); |
| 596 |
Trajectories[nh]=bvus;
|
|
| 594 |
Trajectory[nh]=bvus;
|
|
| 597 | 595 |
|
| 598 | 596 |
} while ((bvus>=rs)&&(bvus<=bvus0)); |
| 599 | 597 |
|
| 598 |
for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {
|
|
| 599 |
Trajectory[i]=0.e0f; |
|
| 600 |
} |
|
| 601 |
|
|
| 600 | 602 |
int imx=(int)(16*bi); |
| 601 | 603 |
|
| 602 | 604 |
for (int i=0;i<imx;i++) |
| 603 | 605 |
{
|
| 604 |
float zp=0,fp=0;
|
|
| 606 |
float zp=0.,fp=0.;
|
|
| 605 | 607 |
float phi=2.*PI/(float)imx*(float)i; |
| 606 | 608 |
float phd=atanp(cos(phi)*sin(tho),cos(tho)); |
| 607 | 609 |
uint yi=(uint)((float)bi*sin(phi)+bmaxi); |
| ... | ... | |
| 618 | 620 |
|
| 619 | 621 |
if (ni<nh) |
| 620 | 622 |
{
|
| 621 |
r=(Trajectories[ni+1]-Trajectories[ni])*(nr-ni*1.)+Trajectories[ni];
|
|
| 623 |
r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
|
|
| 622 | 624 |
} |
| 623 | 625 |
else |
| 624 | 626 |
{
|
| 625 |
r=Trajectories[ni];
|
|
| 627 |
r=Trajectory[ni];
|
|
| 626 | 628 |
} |
| 627 | 629 |
|
| 628 | 630 |
if ((r<=re)&&(r>=ri)) |
| ... | ... | |
| 643 | 645 |
barrier(CLK_GLOBAL_MEM_FENCE); |
| 644 | 646 |
|
| 645 | 647 |
} |
| 648 |
|
|
| 649 |
__kernel void Original(__global float *zImage,__global float *fImage, |
|
| 650 |
uint Size,float Mass,float InternalRadius, |
|
| 651 |
float ExternalRadius,float Angle, |
|
| 652 |
int Line) |
|
| 653 |
{
|
|
| 654 |
// Integer Impact Parameter Size (half of image) |
|
| 655 |
uint bmaxi=(uint)Size; |
|
| 656 |
|
|
| 657 |
float Trajectory[TRACKPOINTS]; |
|
| 658 |
|
|
| 659 |
// Perform trajectory for each pixel |
|
| 660 |
|
|
| 661 |
float m,rs,ri,re,tho; |
|
| 662 |
int raie,q; |
|
| 663 |
|
|
| 664 |
m=Mass; |
|
| 665 |
rs=2.*m; |
|
| 666 |
ri=InternalRadius; |
|
| 667 |
re=ExternalRadius; |
|
| 668 |
tho=Angle; |
|
| 669 |
q=-2; |
|
| 670 |
raie=Line; |
|
| 671 |
|
|
| 672 |
float bmx,db,b,h; |
|
| 673 |
int nh; |
|
| 674 |
|
|
| 675 |
// Autosize for image |
|
| 676 |
bmx=1.25e0f*re; |
|
| 677 |
|
|
| 678 |
// Angular step of integration |
|
| 679 |
h=4.e0f*PI/(float)TRACKPOINTS; |
|
| 680 |
|
|
| 681 |
// Integer Impact Parameter ID |
|
| 682 |
for (int bi=0;bi<bmaxi;bi++) |
|
| 683 |
{
|
|
| 684 |
// impact parameter |
|
| 685 |
b=(float)bi/(float)bmaxi*bmx; |
|
| 686 |
db=bmx/(2.e0f*(float)bmaxi); |
|
| 687 |
|
|
| 688 |
float up,vp,pp,us,vs,ps; |
|
| 689 |
|
|
| 690 |
up=0.; |
|
| 691 |
vp=1.; |
|
| 692 |
|
|
| 693 |
pp=0.; |
|
| 694 |
nh=0; |
|
| 695 |
|
|
| 696 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 697 |
|
|
| 698 |
// b versus us |
|
| 699 |
float bvus=fabs(b/us); |
|
| 700 |
float bvus0=bvus; |
|
| 701 |
Trajectory[nh]=bvus; |
|
| 702 |
|
|
| 703 |
do |
|
| 704 |
{
|
|
| 705 |
nh++; |
|
| 706 |
pp=ps; |
|
| 707 |
up=us; |
|
| 708 |
vp=vs; |
|
| 709 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 710 |
bvus=fabs(b/us); |
|
| 711 |
Trajectory[nh]=bvus; |
|
| 712 |
|
|
| 713 |
} while ((bvus>=rs)&&(bvus<=bvus0)); |
|
| 714 |
|
|
| 715 |
for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {
|
|
| 716 |
Trajectory[i]=0.e0f; |
|
| 717 |
} |
|
| 718 |
|
|
| 719 |
int imx=(int)(16*bi); |
|
| 720 |
|
|
| 721 |
for (int i=0;i<imx;i++) |
|
| 722 |
{
|
|
| 723 |
float zp=0,fp=0; |
|
| 724 |
float phi=2.e0f*PI/(float)imx*(float)i; |
|
| 725 |
float phd=atanp(cos(phi)*sin(tho),cos(tho)); |
|
| 726 |
uint yi=(uint)((float)bi*sin(phi)+bmaxi); |
|
| 727 |
uint xi=(uint)((float)bi*cos(phi)+bmaxi); |
|
| 728 |
|
|
| 729 |
int HalfLap=0,ExitOnImpact=0,ni; |
|
| 730 |
float php,nr,r; |
|
| 731 |
|
|
| 732 |
do |
|
| 733 |
{
|
|
| 734 |
php=phd+(float)HalfLap*PI; |
|
| 735 |
nr=php/h; |
|
| 736 |
ni=(int)nr; |
|
| 737 |
|
|
| 738 |
if (ni<nh) |
|
| 739 |
{
|
|
| 740 |
r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni]; |
|
| 741 |
} |
|
| 742 |
else |
|
| 743 |
{
|
|
| 744 |
r=Trajectory[ni]; |
|
| 745 |
} |
|
| 746 |
|
|
| 747 |
if ((r<=re)&&(r>=ri)) |
|
| 748 |
{
|
|
| 749 |
ExitOnImpact=1; |
|
| 750 |
impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie); |
|
| 751 |
} |
|
| 752 |
|
|
| 753 |
HalfLap++; |
|
| 754 |
|
|
| 755 |
} while ((HalfLap<=2)&&(ExitOnImpact==0)); |
|
| 756 |
|
|
| 757 |
zImage[yi+2*bmaxi*xi]=zp; |
|
| 758 |
fImage[yi+2*bmaxi*xi]=fp; |
|
| 759 |
|
|
| 760 |
} |
|
| 761 |
|
|
| 762 |
} |
|
| 763 |
|
|
| 764 |
barrier(CLK_GLOBAL_MEM_FENCE); |
|
| 765 |
|
|
| 766 |
} |
|
| 646 | 767 |
""" |
| 647 | 768 |
|
| 648 | 769 |
def KernelCodeCuda(): |
| ... | ... | |
| 1145 | 1266 |
|
| 1146 | 1267 |
} |
| 1147 | 1268 |
|
| 1148 |
__global__ void Original(float *zImage,float *fImage,
|
|
| 1149 |
float Mass,float InternalRadius, |
|
| 1150 |
float ExternalRadius,float Angle, |
|
| 1151 |
int Line) |
|
| 1269 |
__global__ void EachCircle(float *zImage,float *fImage,
|
|
| 1270 |
float Mass,float InternalRadius,
|
|
| 1271 |
float ExternalRadius,float Angle,
|
|
| 1272 |
int Line)
|
|
| 1152 | 1273 |
{
|
| 1153 | 1274 |
// Integer Impact Parameter ID |
| 1154 | 1275 |
int bi=blockIdx.x*blockDim.x+threadIdx.x; |
| ... | ... | |
| 1156 | 1277 |
// Integer Impact Parameter Size (half of image) |
| 1157 | 1278 |
int bmaxi=gridDim.x*blockDim.x; |
| 1158 | 1279 |
|
| 1159 |
float Trajectories[2048];
|
|
| 1280 |
float Trajectory[2048];
|
|
| 1160 | 1281 |
|
| 1161 | 1282 |
// Perform trajectory for each pixel |
| 1162 | 1283 |
|
| ... | ... | |
| 1197 | 1318 |
// b versus us |
| 1198 | 1319 |
float bvus=fabs(b/us); |
| 1199 | 1320 |
float bvus0=bvus; |
| 1200 |
Trajectories[nh]=bvus;
|
|
| 1321 |
Trajectory[nh]=bvus;
|
|
| 1201 | 1322 |
|
| 1202 | 1323 |
do |
| 1203 | 1324 |
{
|
| ... | ... | |
| 1207 | 1328 |
vp=vs; |
| 1208 | 1329 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
| 1209 | 1330 |
bvus=fabs(b/us); |
| 1210 |
Trajectories[nh]=bvus;
|
|
| 1331 |
Trajectory[nh]=bvus;
|
|
| 1211 | 1332 |
|
| 1212 | 1333 |
} while ((bvus>=rs)&&(bvus<=bvus0)); |
| 1213 | 1334 |
|
| ... | ... | |
| 1232 | 1353 |
|
| 1233 | 1354 |
if (ni<nh) |
| 1234 | 1355 |
{
|
| 1235 |
r=(Trajectories[ni+1]-Trajectories[ni])*(nr-ni*1.)+Trajectories[ni];
|
|
| 1356 |
r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
|
|
| 1236 | 1357 |
} |
| 1237 | 1358 |
else |
| 1238 | 1359 |
{
|
| 1239 |
r=Trajectories[ni];
|
|
| 1360 |
r=Trajectory[ni];
|
|
| 1240 | 1361 |
} |
| 1241 | 1362 |
|
| 1242 | 1363 |
if ((r<=re)&&(r>=ri)) |
| ... | ... | |
| 1257 | 1378 |
} |
| 1258 | 1379 |
|
| 1259 | 1380 |
} |
| 1381 |
|
|
| 1382 |
__global__ void Original(float *zImage,float *fImage, |
|
| 1383 |
uint Size,float Mass,float InternalRadius, |
|
| 1384 |
float ExternalRadius,float Angle, |
|
| 1385 |
int Line) |
|
| 1386 |
{
|
|
| 1387 |
// Integer Impact Parameter Size (half of image) |
|
| 1388 |
uint bmaxi=(uint)Size; |
|
| 1389 |
|
|
| 1390 |
float Trajectory[TRACKPOINTS]; |
|
| 1391 |
|
|
| 1392 |
// Perform trajectory for each pixel |
|
| 1393 |
|
|
| 1394 |
float m,rs,ri,re,tho; |
|
| 1395 |
int raie,q; |
|
| 1396 |
|
|
| 1397 |
m=Mass; |
|
| 1398 |
rs=2.*m; |
|
| 1399 |
ri=InternalRadius; |
|
| 1400 |
re=ExternalRadius; |
|
| 1401 |
tho=Angle; |
|
| 1402 |
q=-2; |
|
| 1403 |
raie=Line; |
|
| 1404 |
|
|
| 1405 |
float bmx,db,b,h; |
|
| 1406 |
int nh; |
|
| 1407 |
|
|
| 1408 |
// Autosize for image |
|
| 1409 |
bmx=1.25e0f*re; |
|
| 1410 |
|
|
| 1411 |
// Angular step of integration |
|
| 1412 |
h=4.e0f*PI/(float)TRACKPOINTS; |
|
| 1413 |
|
|
| 1414 |
// Integer Impact Parameter ID |
|
| 1415 |
for (int bi=0;bi<bmaxi;bi++) |
|
| 1416 |
{
|
|
| 1417 |
// impact parameter |
|
| 1418 |
b=(float)bi/(float)bmaxi*bmx; |
|
| 1419 |
db=bmx/(2.e0f*(float)bmaxi); |
|
| 1420 |
|
|
| 1421 |
float up,vp,pp,us,vs,ps; |
|
| 1422 |
|
|
| 1423 |
up=0.; |
|
| 1424 |
vp=1.; |
|
| 1425 |
|
|
| 1426 |
pp=0.; |
|
| 1427 |
nh=0; |
|
| 1428 |
|
|
| 1429 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 1430 |
|
|
| 1431 |
// b versus us |
|
| 1432 |
float bvus=fabs(b/us); |
|
| 1433 |
float bvus0=bvus; |
|
| 1434 |
Trajectory[nh]=bvus; |
|
| 1435 |
|
|
| 1436 |
do |
|
| 1437 |
{
|
|
| 1438 |
nh++; |
|
| 1439 |
pp=ps; |
|
| 1440 |
up=us; |
|
| 1441 |
vp=vs; |
|
| 1442 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 1443 |
bvus=fabs(b/us); |
|
| 1444 |
Trajectory[nh]=bvus; |
|
| 1445 |
|
|
| 1446 |
} while ((bvus>=rs)&&(bvus<=bvus0)); |
|
| 1447 |
|
|
| 1448 |
for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {
|
|
| 1449 |
Trajectory[i]=0.e0f; |
|
| 1450 |
} |
|
| 1451 |
|
|
| 1452 |
int imx=(int)(16*bi); |
|
| 1453 |
|
|
| 1454 |
for (int i=0;i<imx;i++) |
|
| 1455 |
{
|
|
| 1456 |
float zp=0,fp=0; |
|
| 1457 |
float phi=2.e0f*PI/(float)imx*(float)i; |
|
| 1458 |
float phd=atanp(cos(phi)*sin(tho),cos(tho)); |
|
| 1459 |
uint yi=(uint)((float)bi*sin(phi)+bmaxi); |
|
| 1460 |
uint xi=(uint)((float)bi*cos(phi)+bmaxi); |
|
| 1461 |
|
|
| 1462 |
int HalfLap=0,ExitOnImpact=0,ni; |
|
| 1463 |
float php,nr,r; |
|
| 1464 |
|
|
| 1465 |
do |
|
| 1466 |
{
|
|
| 1467 |
php=phd+(float)HalfLap*PI; |
|
| 1468 |
nr=php/h; |
|
| 1469 |
ni=(int)nr; |
|
| 1470 |
|
|
| 1471 |
if (ni<nh) |
|
| 1472 |
{
|
|
| 1473 |
r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni]; |
|
| 1474 |
} |
|
| 1475 |
else |
|
| 1476 |
{
|
|
| 1477 |
r=Trajectory[ni]; |
|
| 1478 |
} |
|
| 1479 |
|
|
| 1480 |
if ((r<=re)&&(r>=ri)) |
|
| 1481 |
{
|
|
| 1482 |
ExitOnImpact=1; |
|
| 1483 |
impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie); |
|
| 1484 |
} |
|
| 1485 |
|
|
| 1486 |
HalfLap++; |
|
| 1487 |
|
|
| 1488 |
} while ((HalfLap<=2)&&(ExitOnImpact==0)); |
|
| 1489 |
|
|
| 1490 |
zImage[yi+2*bmaxi*xi]=zp; |
|
| 1491 |
fImage[yi+2*bmaxi*xi]=fp; |
|
| 1492 |
|
|
| 1493 |
} |
|
| 1494 |
|
|
| 1495 |
} |
|
| 1496 |
|
|
| 1497 |
} |
|
| 1260 | 1498 |
""" |
| 1261 | 1499 |
return(BlobCUDA) |
| 1262 | 1500 |
|
| ... | ... | |
| 1356 | 1594 |
numpy.int32(Line)) |
| 1357 | 1595 |
CLLaunch.wait() |
| 1358 | 1596 |
elif Method=='Original': |
| 1359 |
CLLaunch=BlackHoleCL.Original(queue,(zImage.shape[0]/2,),
|
|
| 1597 |
CLLaunch=BlackHoleCL.Original(queue,(1,),
|
|
| 1360 | 1598 |
None,zImageCL,fImageCL, |
| 1599 |
numpy.uint32(zImage.shape[0]/2), |
|
| 1361 | 1600 |
numpy.float32(Mass), |
| 1362 | 1601 |
numpy.float32(InternalRadius), |
| 1363 | 1602 |
numpy.float32(ExternalRadius), |
| 1364 | 1603 |
numpy.float32(Angle), |
| 1365 | 1604 |
numpy.int32(Line)) |
| 1366 | 1605 |
CLLaunch.wait() |
| 1606 |
elif Method=='EachCircle': |
|
| 1607 |
CLLaunch=BlackHoleCL.EachCircle(queue,(zImage.shape[0]/2,), |
|
| 1608 |
None,zImageCL,fImageCL, |
|
| 1609 |
numpy.float32(Mass), |
|
| 1610 |
numpy.float32(InternalRadius), |
|
| 1611 |
numpy.float32(ExternalRadius), |
|
| 1612 |
numpy.float32(Angle), |
|
| 1613 |
numpy.int32(Line)) |
|
| 1614 |
CLLaunch.wait() |
|
| 1367 | 1615 |
elif Method=='TrajectoCircle': |
| 1368 | 1616 |
CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],), |
| 1369 | 1617 |
None,TrajectoriesCL,IdLastCL, |
| ... | ... | |
| 1480 | 1728 |
|
| 1481 | 1729 |
EachPixelCU=mod.get_function("EachPixel")
|
| 1482 | 1730 |
OriginalCU=mod.get_function("Original")
|
| 1731 |
EachCircleCU=mod.get_function("EachCircle")
|
|
| 1483 | 1732 |
TrajectoryCU=mod.get_function("Trajectory")
|
| 1484 | 1733 |
PixelCU=mod.get_function("Pixel")
|
| 1485 | 1734 |
CircleCU=mod.get_function("Circle")
|
| ... | ... | |
| 1504 | 1753 |
numpy.int32(Line), |
| 1505 | 1754 |
grid=(zImage.shape[0]/32,zImage.shape[1]/32), |
| 1506 | 1755 |
block=(32,32,1)) |
| 1756 |
elif Method=='EachCircle': |
|
| 1757 |
EachCircleCU(zImageCU,fImageCU, |
|
| 1758 |
numpy.float32(Mass), |
|
| 1759 |
numpy.float32(InternalRadius), |
|
| 1760 |
numpy.float32(ExternalRadius), |
|
| 1761 |
numpy.float32(Angle), |
|
| 1762 |
numpy.int32(Line), |
|
| 1763 |
grid=(zImage.shape[0]/32/2,1), |
|
| 1764 |
block=(32,1,1)) |
|
| 1507 | 1765 |
elif Method=='Original': |
| 1508 | 1766 |
OriginalCU(zImageCU,fImageCU, |
| 1767 |
numpy.uint32(zImage.shape[0]/2), |
|
| 1509 | 1768 |
numpy.float32(Mass), |
| 1510 | 1769 |
numpy.float32(InternalRadius), |
| 1511 | 1770 |
numpy.float32(ExternalRadius), |
| 1512 | 1771 |
numpy.float32(Angle), |
| 1513 | 1772 |
numpy.int32(Line), |
| 1514 |
grid=(zImage.shape[0]/32/2,1),
|
|
| 1515 |
block=(32,1,1))
|
|
| 1773 |
grid=(1,1),
|
|
| 1774 |
block=(1,1,1))
|
|
| 1516 | 1775 |
elif Method=='TrajectoCircle': |
| 1517 | 1776 |
TrajectoryCU(TrajectoriesCU,IdLastCU, |
| 1518 | 1777 |
numpy.float32(Mass), |
| ... | ... | |
| 1601 | 1860 |
# No output as image |
| 1602 | 1861 |
NoImage = False |
| 1603 | 1862 |
|
| 1604 |
HowToUse='%s -h [Help] -b [BlackBodyEmission] -n [NoImage] -p <Einstein/Newton> -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -a <AngleAboveDisc> -d <DeviceId> -c <Greyscale/Red2Yellow> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel/Original> -v <FP32/FP64>' |
|
| 1863 |
HowToUse='%s -h [Help] -b [BlackBodyEmission] -n [NoImage] -p <Einstein/Newton> -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -a <AngleAboveDisc> -d <DeviceId> -c <Greyscale/Red2Yellow> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel/EachCircle/Original> -v <FP32/FP64>'
|
|
| 1605 | 1864 |
|
| 1606 | 1865 |
try: |
| 1607 | 1866 |
opts, args = getopt.getopt(sys.argv[1:],"hbns:m:i:x:a:d:g:v:t:c:p:",["blackbody","noimage","camera","size=","mass=","internal=","external=","angle=","device=","gpustyle=","variabletype=","method=","colors=","physics="]) |
Formats disponibles : Unified diff