Révision 211 TrouNoir/TrouNoir.py
| TrouNoir.py (revision 211) | ||
|---|---|---|
| 26 | 26 |
|
| 27 | 27 |
import pyopencl as cl |
| 28 | 28 |
import numpy |
| 29 |
from PIL import Image |
|
| 30 | 29 |
import time,string |
| 31 | 30 |
from numpy.random import randint as nprnd |
| 32 | 31 |
import sys |
| 33 | 32 |
import getopt |
| 34 | 33 |
import matplotlib.pyplot as plt |
| 34 |
from socket import gethostname |
|
| 35 | 35 |
|
| 36 |
def DictionariesAPI(): |
|
| 37 |
PhysicsList={'Einstein':0,'Newton':1}
|
|
| 38 |
return(PhysicsList) |
|
| 36 | 39 |
|
| 40 |
BlobOpenCL= """ |
|
| 37 | 41 |
|
| 38 |
|
|
| 39 |
|
|
| 40 |
|
|
| 41 |
|
|
| 42 |
|
|
| 43 |
KERNEL_CODE=string.Template("""
|
|
| 44 |
|
|
| 45 | 42 |
#define PI (float)3.14159265359 |
| 46 | 43 |
#define nbr 256 |
| 47 | 44 |
|
| 45 |
#define EINSTEIN 0 |
|
| 46 |
#define NEWTON 1 |
|
| 47 |
|
|
| 48 | 48 |
float atanp(float x,float y) |
| 49 | 49 |
{
|
| 50 | 50 |
float angle; |
| ... | ... | |
| 64 | 64 |
return v; |
| 65 | 65 |
} |
| 66 | 66 |
|
| 67 |
#if PHYSICS == NEWTON |
|
| 67 | 68 |
float g(float u,float m,float b) |
| 68 | 69 |
{
|
| 70 |
return (-u); |
|
| 71 |
} |
|
| 72 |
#else |
|
| 73 |
float g(float u,float m,float b) |
|
| 74 |
{
|
|
| 69 | 75 |
return (3.e0f*m/b*pow(u,2)-u); |
| 70 | 76 |
} |
| 77 |
#endif |
|
| 71 | 78 |
|
| 72 |
|
|
| 73 | 79 |
void calcul(float *us,float *vs,float up,float vp, |
| 74 | 80 |
float h,float m,float b) |
| 75 | 81 |
{
|
| ... | ... | |
| 156 | 162 |
if ((posfreq>0)&&(posfreq<nbr)) |
| 157 | 163 |
{
|
| 158 | 164 |
// Initial version |
| 159 |
flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.); |
|
| 165 |
// flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.);
|
|
| 160 | 166 |
// Version with log used |
| 161 |
flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.); |
|
| 162 |
flux_int*=pow(rf,3)*b*db*h; |
|
| 167 |
//flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.); |
|
| 168 |
// flux_int*=pow(rf,3)*b*db*h; |
|
| 169 |
//flux_int*=exp(3.*log(rf))*b*db*h; |
|
| 170 |
flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.)*exp(3.*log(rf))*b*db*h; |
|
| 171 |
|
|
| 163 | 172 |
flx+=flux_int; |
| 164 | 173 |
} |
| 165 | 174 |
} |
| ... | ... | |
| 522 | 531 |
barrier(CLK_GLOBAL_MEM_FENCE); |
| 523 | 532 |
|
| 524 | 533 |
} |
| 525 |
""")
|
|
| 534 |
""" |
|
| 526 | 535 |
|
| 527 |
def ImageOutput(sigma,prefix): |
|
| 528 |
Max=sigma.max() |
|
| 529 |
Min=sigma.min() |
|
| 536 |
# def ImageOutput(sigma,prefix): |
|
| 537 |
# from PIL import Image |
|
| 538 |
# Max=sigma.max() |
|
| 539 |
# Min=sigma.min() |
|
| 530 | 540 |
|
| 531 |
# Normalize value as 8bits Integer |
|
| 532 |
SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
|
|
| 533 |
image = Image.fromarray(SigmaInt) |
|
| 534 |
image.save("%s.jpg" % prefix)
|
|
| 535 |
|
|
| 541 |
# # Normalize value as 8bits Integer |
|
| 542 |
# SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8')
|
|
| 543 |
# image = Image.fromarray(SigmaInt) |
|
| 544 |
# image.save("%s.jpg" % prefix)
|
|
| 545 |
|
|
| 546 |
def ImageOutput(sigma,prefix,Colors): |
|
| 547 |
import matplotlib.pyplot as plt |
|
| 548 |
if Colors == 'Red2Yellow': |
|
| 549 |
plt.imsave("%s.png" % prefix, sigma, cmap='afmhot')
|
|
| 550 |
else: |
|
| 551 |
plt.imsave("%s.png" % prefix, sigma, cmap='Greys_r')
|
|
| 552 |
|
|
| 536 | 553 |
def BlackHoleCL(zImage,fImage,InputCL): |
| 537 | 554 |
|
| 538 | 555 |
kernel_params = {}
|
| ... | ... | |
| 549 | 566 |
Angle=InputCL['Angle'] |
| 550 | 567 |
Method=InputCL['Method'] |
| 551 | 568 |
TrackPoints=InputCL['TrackPoints'] |
| 569 |
Physics=InputCL['Physics'] |
|
| 552 | 570 |
|
| 571 |
PhysicsList=DictionariesAPI() |
|
| 572 |
|
|
| 553 | 573 |
if InputCL['BlackBody']: |
| 554 | 574 |
# Spectrum is Black Body one |
| 555 | 575 |
Line=0 |
| ... | ... | |
| 580 | 600 |
queue = cl.CommandQueue(ctx, |
| 581 | 601 |
properties=cl.command_queue_properties.PROFILING_ENABLE) |
| 582 | 602 |
|
| 583 |
BlackHoleCL = cl.Program(ctx,KERNEL_CODE.substitute(kernel_params)).build() |
|
| 584 | 603 |
|
| 604 |
# BlackHoleCL = cl.Program(ctx,KERNEL_CODE.substitute(kernel_params)).build() |
|
| 605 |
|
|
| 606 |
BuildOptions="-cl-mad-enable -DPHYSICS=%i " % (PhysicsList[Physics]) |
|
| 607 |
|
|
| 608 |
BlackHoleCL = cl.Program(ctx,BlobOpenCL).build(options = BuildOptions) |
|
| 609 |
|
|
| 585 | 610 |
# Je recupere les flag possibles pour les buffers |
| 586 | 611 |
mf = cl.mem_flags |
| 587 | 612 |
|
| 588 |
print(zImage) |
|
| 589 |
|
|
| 590 | 613 |
TrajectoriesCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=Trajectories) |
| 591 | 614 |
zImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=zImage) |
| 592 | 615 |
fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage) |
| ... | ... | |
| 594 | 617 |
|
| 595 | 618 |
start_time=time.time() |
| 596 | 619 |
|
| 597 |
print(Trajectories.shape) |
|
| 598 | 620 |
if Method=='EachPixel': |
| 599 | 621 |
CLLaunch=BlackHoleCL.EachPixel(queue,(zImage.shape[0],zImage.shape[1]), |
| 600 | 622 |
None,zImageCL,fImageCL, |
| ... | ... | |
| 648 | 670 |
CLLaunch.wait() |
| 649 | 671 |
|
| 650 | 672 |
elapsed = time.time()-start_time |
| 651 |
print("Elapsed %f: " % elapsed)
|
|
| 673 |
print("\nElapsed Time : %f" % elapsed)
|
|
| 652 | 674 |
|
| 653 | 675 |
cl.enqueue_copy(queue,zImage,zImageCL).wait() |
| 654 | 676 |
cl.enqueue_copy(queue,fImage,fImageCL).wait() |
| 655 | 677 |
cl.enqueue_copy(queue,Trajectories,TrajectoriesCL).wait() |
| 656 | 678 |
cl.enqueue_copy(queue,IdLast,IdLastCL).wait() |
| 657 |
|
|
| 658 |
print(zImage.max(),numpy.where(zImage[:,:]==zImage.max())) |
|
| 659 |
print(fImage.max(),numpy.where(fImage[:,:]==fImage.max())) |
|
| 679 |
|
|
| 680 |
zMaxPosition=numpy.where(zImage[:,:]==zImage.max()) |
|
| 681 |
fMaxPosition=numpy.where(fImage[:,:]==fImage.max()) |
|
| 682 |
print("Z max @(%i,%i) : %f" % (zMaxPosition[0],zMaxPosition[1],zImage.max()))
|
|
| 683 |
print("Flux max @(%i,%i) : %f\n" % (fMaxPosition[0],fMaxPosition[1],fImage.max()))
|
|
| 660 | 684 |
zImageCL.release() |
| 661 | 685 |
fImageCL.release() |
| 662 | 686 |
|
| ... | ... | |
| 686 | 710 |
q=-2 |
| 687 | 711 |
# Method of resolution |
| 688 | 712 |
Method='TrajectoPixel' |
| 689 |
|
|
| 690 |
HowToUse='%s -h [Help] -b [BlackBodyEmission] -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -a <AngleAboveDisc> -d <DeviceId> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel> -v <FP32/FP64>' |
|
| 713 |
# Colors for output image |
|
| 714 |
Colors='Greyscale' |
|
| 715 |
# Physics |
|
| 716 |
Physics='Einstein' |
|
| 717 |
# No output as image |
|
| 718 |
NoImage = False |
|
| 719 |
|
|
| 720 |
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> -v <FP32/FP64>' |
|
| 691 | 721 |
|
| 692 | 722 |
try: |
| 693 |
opts, args = getopt.getopt(sys.argv[1:],"hbs:m:i:x:a:d:g:v:t:",["blackbody","camera","size=","mass=","internal=","external=","angle=","device=","gpustyle=","variabletype=","method="])
|
|
| 723 |
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="])
|
|
| 694 | 724 |
except getopt.GetoptError: |
| 695 | 725 |
print(HowToUse % sys.argv[0]) |
| 696 | 726 |
sys.exit(2) |
| ... | ... | |
| 751 | 781 |
Angle = numpy.pi/180.*(90.-float(arg)) |
| 752 | 782 |
elif opt in ("-b", "--blackbody"):
|
| 753 | 783 |
BlackBody = True |
| 784 |
elif opt in ("-n", "--noimage"):
|
|
| 785 |
NoImage = True |
|
| 754 | 786 |
elif opt in ("-t", "--method"):
|
| 755 | 787 |
Method = arg |
| 788 |
elif opt in ("-c", "--colors"):
|
|
| 789 |
Colors = arg |
|
| 790 |
elif opt in ("-p", "--physics"):
|
|
| 791 |
Physics = arg |
|
| 756 | 792 |
|
| 757 | 793 |
print("Device Identification selected : %s" % Device)
|
| 758 | 794 |
print("GpuStyle used : %s" % GpuStyle)
|
| ... | ... | |
| 764 | 800 |
print("Angle with normal of (in radians) : %f" % Angle)
|
| 765 | 801 |
print("Black Body Disc Emission (monochromatic instead) : %s" % BlackBody)
|
| 766 | 802 |
print("Method of resolution : %s" % Method)
|
| 803 |
print("Colors for output images : %s" % Colors)
|
|
| 804 |
print("Physics used for Trajectories : %s" % Physics)
|
|
| 767 | 805 |
|
| 768 | 806 |
if GpuStyle=='CUDA': |
| 769 | 807 |
print("\nSelection of CUDA device")
|
| ... | ... | |
| 819 | 857 |
InputCL['BlackBody']=BlackBody |
| 820 | 858 |
InputCL['Method']=Method |
| 821 | 859 |
InputCL['TrackPoints']=TrackPoints |
| 860 |
InputCL['Physics']=Physics |
|
| 822 | 861 |
|
| 823 | 862 |
duration=BlackHoleCL(zImage,fImage,InputCL) |
| 824 | 863 |
|
| 825 |
ImageOutput(zImage,"TrouNoirZ_%s" % Method) |
|
| 826 |
ImageOutput(fImage,"TrouNoirF_%s" % Method) |
|
| 864 |
Hostname=gethostname() |
|
| 865 |
Date=time.strftime("%Y%m%d_%H%M%S")
|
|
| 866 |
ImageInfo="%s_Device%i_%s_%s" % (Method,Device,Hostname,Date) |
|
| 867 |
|
|
| 868 |
|
|
| 869 |
if not NoImage: |
|
| 870 |
ImageOutput(zImage,"TrouNoirZ_%s" % ImageInfo,Colors) |
|
| 871 |
ImageOutput(fImage,"TrouNoirF_%s" % ImageInfo,Colors) |
|
Formats disponibles : Unified diff