#!/usr/bin/env python # # TrouNoir model using PyOpenCL # # CC BY-NC-SA 2019 : # # Thanks to Andreas Klockner for PyOpenCL: # http://mathema.tician.de/software/pyopencl # # Original code programmed in Fortran 77 in mars 1994 # for Practical Work of Numerical Simulation # DEA (old Master2) in astrophysics and spatial techniques in Meudon # by Herve Aussel & Emmanuel Quemener # # Conversion in C done by Emmanuel Quemener in august 1997 # GPUfication in OpenCL under Python in july 2019 # # Thanks to : # # - Herve Aussel for his part of code of black body spectrum # - Didier Pelat for his help to perform this work # - Jean-Pierre Luminet for his article published in 1979 # - Numerical Recipies for Runge Kutta recipies # - Luc Blanchet for his disponibility about my questions in General Relativity # - Pierre Lena for his passion about science and vulgarisation import pyopencl as cl import numpy from PIL import Image import time,string from numpy.random import randint as nprnd import sys import getopt import matplotlib.pyplot as plt KERNEL_CODE=string.Template(""" #define PI (float)3.14159265359 #define nbr 200 float atanp(float x,float y) { float angle; angle=atan2(y,x); if (angle<0.e0f) { angle+=(float)2.e0f*PI; } return angle; } float f(float v) { return v; } float g(float u,float m,float b) { return (3.e0f*m/b*pow(u,2)-u); } void calcul(float *us,float *vs,float up,float vp, float h,float m,float b) { float c0,c1,c2,c3,d0,d1,d2,d3; c0=h*f(vp); c1=h*f(vp+c0/2.); c2=h*f(vp+c1/2.); c3=h*f(vp+c2); d0=h*g(up,m,b); d1=h*g(up+d0/2.,m,b); d2=h*g(up+d1/2.,m,b); d3=h*g(up+d2,m,b); *us=up+(c0+2.*c1+2.*c2+c3)/6.; *vs=vp+(d0+2.*d1+2.*d2+d3)/6.; } void rungekutta(float *ps,float *us,float *vs, float pp,float up,float vp, float h,float m,float b) { calcul(us,vs,up,vp,h,m,b); *ps=pp+h; } float decalage_spectral(float r,float b,float phi, float tho,float m) { return (sqrt(1-3*m/r)/(1+sqrt(m/pow(r,3))*b*sin(tho)*sin(phi))); } float spectre(float rf,int q,float b,float db, float h,float r,float m,float bss) { float flx; flx=pow(r/m,q)*pow(rf,4)*b; return(flx); } float spectre_cn(float rf,float b,float db,float h,float r,float m,float bss) { double flx; double nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k,psc2,psk; int fi,posfreq; planck=6.62e-34; k=1.38e-23; psk=5.38e-11; psc2=7.35e-51; temp=3.e7; m_point=1.e14; v=1.-3./r; qu=1./sqrt(1.-3./r)/sqrt(r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))*(sqrt(6.)-sqrt(3.))/(sqrt(6.)+sqrt(3.)))); temp_em=temp*sqrt(m)*exp(0.25*log(m_point))*exp(-0.75*log(r))*exp(-0.125*log(v))*exp(0.25*log(qu)); flux_int=0; flx=0; for (fi=0;fi0)&&(posfreqfmod(phd,PI)))&&(rp[nh]>ri)&&(rp[nh]=rs)&&(rp[nh]<=rp[0])&&(ExitOnImpact==0)); if (ExitOnImpact==1) { impact(phi,rp[nh-1],b,tho,m,&zp,&fp,q,db,h,raie); } else { zp=0.; fp=0.; } barrier(CLK_GLOBAL_MEM_FENCE); zImage[yi+sizex*xi]=(float)zp; fImage[yi+sizex*xi]=(float)fp; } __kernel void Pixel(__global float *zImage,__global float *fImage, __global float *Trajectories,__global int *IdLast, uint ImpactParameter,uint TrackPoints, float Mass,float InternalRadius, float ExternalRadius,float Angle, float ObserverDistance, int BlackBody,int AngularCamera) { uint xi=(uint)get_global_id(0); uint yi=(uint)get_global_id(1); uint sizex=(uint)get_global_size(0); uint sizey=(uint)get_global_size(1); // Perform trajectory for each pixel float m,rs,ri,re,tho,ro; int q,dist,raie; m=Mass; rs=2.*m; ri=InternalRadius; re=ExternalRadius; ro=ObserverDistance; tho=Angle; q=-2; dist=AngularCamera; raie=BlackBody; float d,bmx,db,b,h; float phi,thi,phd,php,nr,r; int nh; float zp=0,fp=0; // Autosize for image, 25% greater than externa radius bmx=1.25*re; // Angular step of integration h=4.e0f*PI/(float)TrackPoints; // Step of Impact Parameter db=bmx/(float)ImpactParameter; // set origin as center of image float x=(float)xi-(float)(sizex/2)+(float)5e-1f; float y=(float)yi-(float)(sizey/2)+(float)5e-1f; // angle extracted from cylindric symmetry phi=atanp(x,y); phd=atanp(cos(phi)*sin(tho),cos(tho)); // Real Impact Parameter b=sqrt(x*x+y*y)*bmx/(float)ImpactParameter; // Integer Impact Parameter uint bi=(uint)sqrt(x*x+y*y); int HalfLap=0,ExitOnImpact=0,ni; if (bi=ri)) { ExitOnImpact=1; impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie); } HalfLap++; } while ((HalfLap<=2)&&(ExitOnImpact==0)); } barrier(CLK_GLOBAL_MEM_FENCE); zImage[yi+sizex*xi]=zp; fImage[yi+sizex*xi]=fp; } __kernel void Circle(__global float *Trajectories,__global int *IdLast, __global float *zImage,__global float *fImage, int TrackPoints, float Mass,float InternalRadius, float ExternalRadius,float Angle, float ObserverDistance, int BlackBody,int AngularCamera) { // Integer Impact Parameter ID int bi=get_global_id(0); // Integer points on circle int i=get_global_id(1); // Integer Impact Parameter Size (half of image) int bmaxi=get_global_size(0); // Integer Points on circle int imx=get_global_size(1); // Perform trajectory for each pixel float m,rs,ri,re,tho,ro; int q,dist,raie; m=Mass; rs=2.*m; ri=InternalRadius; re=ExternalRadius; ro=ObserverDistance; tho=Angle; q=-2; dist=AngularCamera; raie=BlackBody; float bmx,db,b,h; float phi,thi,phd; int nh; float zp=0,fp=0; // Autosize for image bmx=1.25*re; // Angular step of integration h=4.e0f*PI/(float)TrackPoints; // impact parameter b=(float)bi/(float)bmaxi*bmx; db=bmx/(float)bmaxi; phi=2.*PI/(float)imx*(float)i; phd=atanp(cos(phi)*sin(tho),cos(tho)); int yi=(int)((float)bi*sin(phi))+bmaxi; int xi=(int)((float)bi*cos(phi))+bmaxi; int HalfLap=0,ExitOnImpact=0,ni; float php,nr,r; do { php=phd+(float)HalfLap*PI; nr=php/h; ni=(int)nr; if (ni=ri)) { ExitOnImpact=1; impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie); } HalfLap++; } while ((HalfLap<=2)&&(ExitOnImpact==0)); zImage[yi+2*bmaxi*xi]=zp; fImage[yi+2*bmaxi*xi]=fp; barrier(CLK_GLOBAL_MEM_FENCE); } __kernel void Trajectory(__global float *Trajectories, __global int *IdLast,int TrackPoints, float Mass,float InternalRadius, float ExternalRadius,float Angle, float ObserverDistance, int BlackBody,int AngularCamera) { // Integer Impact Parameter ID int bi=get_global_id(0); // Integer Impact Parameter Size (half of image) int bmaxi=get_global_size(0); // Perform trajectory for each pixel float m,rs,ri,re,tho,ro; int dist,raie,q; m=Mass; rs=2.*m; ri=InternalRadius; re=ExternalRadius; ro=ObserverDistance; tho=Angle; q=-2; dist=AngularCamera; raie=BlackBody; float d,bmx,db,b,h; float phi,thi,phd,php,nr,r; int nh; float zp,fp; // Autosize for image bmx=1.25*re; // Angular step of integration h=4.e0f*PI/(float)TrackPoints; // impact parameter b=(float)bi/(float)bmaxi*bmx; float up,vp,pp,us,vs,ps; if (dist==1) { up=0.; vp=1.; } else { thi=asin(b/ro); up=sin(thi); vp=cos(thi); } pp=0.; nh=0; rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); // b versus us float bvus=fabs(b/us); float bvus0=bvus; Trajectories[bi*TrackPoints+nh]=bvus; do { nh++; pp=ps; up=us; vp=vs; rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); bvus=fabs(b/us); Trajectories[bi*TrackPoints+nh]=bvus; } while ((bvus>=rs)&&(bvus<=bvus0)); IdLast[bi]=nh; barrier(CLK_GLOBAL_MEM_FENCE); } """) def ImageOutput(sigma,prefix): Max=sigma.max() Min=sigma.min() # Normalize value as 8bits Integer SigmaInt=(255*(sigma-Min)/(Max-Min)).astype('uint8') image = Image.fromarray(SigmaInt) image.save("%s.jpg" % prefix) def BlackHoleCL(zImage,fImage,InputCL): kernel_params = {} print(InputCL) Device=InputCL['Device'] GpuStyle=InputCL['GpuStyle'] VariableType=InputCL['VariableType'] Size=InputCL['Size'] Mass=InputCL['Mass'] InternalRadius=InputCL['InternalRadius'] ExternalRadius=InputCL['ExternalRadius'] Angle=InputCL['Angle'] ObserverDistance=InputCL['ObserverDistance'] Method=InputCL['Method'] TrackPoints=InputCL['TrackPoints'] if InputCL['BlackBody']: BlackBody=0 else: BlackBody=1 if InputCL['AngularCamera']: AngularCamera=1 else: AngularCamera=0 Trajectories=numpy.zeros((InputCL['Size']/2,InputCL['TrackPoints']), dtype=numpy.float32) IdLast=numpy.zeros(InputCL['Size']/2,dtype=numpy.int32) # Je detecte un peripherique GPU dans la liste des peripheriques Id=0 HasXPU=False for platform in cl.get_platforms(): for device in platform.get_devices(): if Id==Device: XPU=device print "CPU/GPU selected: ",device.name.lstrip() HasXPU=True Id+=1 if HasXPU==False: print "No XPU #%i found in all of %i devices, sorry..." % (Device,Id-1) sys.exit() ctx = cl.Context([XPU]) queue = cl.CommandQueue(ctx, properties=cl.command_queue_properties.PROFILING_ENABLE) BlackHoleCL = cl.Program(ctx,KERNEL_CODE.substitute(kernel_params)).build() # Je recupere les flag possibles pour les buffers mf = cl.mem_flags print(zImage) TrajectoriesCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=Trajectories) zImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=zImage) fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage) fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage) IdLastCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=IdLast) start_time=time.time() print(Trajectories.shape) if Method=='EachPixel': CLLaunch=BlackHoleCL.EachPixel(queue,(zImage.shape[0],zImage.shape[1]), None,zImageCL,fImageCL, numpy.float32(Mass), numpy.float32(InternalRadius), numpy.float32(ExternalRadius), numpy.float32(Angle), numpy.float32(ObserverDistance), numpy.int32(BlackBody), numpy.int32(AngularCamera)) CLLaunch.wait() elif Method=='TrajectoCircle': CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],), None,TrajectoriesCL,IdLastCL, numpy.uint32(Trajectories.shape[1]), numpy.float32(Mass), numpy.float32(InternalRadius), numpy.float32(ExternalRadius), numpy.float32(Angle), numpy.float32(ObserverDistance), numpy.int32(BlackBody), numpy.int32(AngularCamera)) CLLaunch=BlackHoleCL.Circle(queue,(Trajectories.shape[0], zImage.shape[0]*4),None, TrajectoriesCL,IdLastCL, zImageCL,fImageCL, numpy.uint32(Trajectories.shape[1]), numpy.float32(Mass), numpy.float32(InternalRadius), numpy.float32(ExternalRadius), numpy.float32(Angle), numpy.float32(ObserverDistance), numpy.int32(BlackBody), numpy.int32(AngularCamera)) CLLaunch.wait() else: CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],), None,TrajectoriesCL,IdLastCL, numpy.uint32(Trajectories.shape[1]), numpy.float32(Mass), numpy.float32(InternalRadius), numpy.float32(ExternalRadius), numpy.float32(Angle), numpy.float32(ObserverDistance), numpy.int32(BlackBody), numpy.int32(AngularCamera)) CLLaunch=BlackHoleCL.Pixel(queue,(zImage.shape[0], zImage.shape[1]),None, zImageCL,fImageCL,TrajectoriesCL,IdLastCL, numpy.uint32(Trajectories.shape[0]), numpy.uint32(Trajectories.shape[1]), numpy.float32(Mass), numpy.float32(InternalRadius), numpy.float32(ExternalRadius), numpy.float32(Angle), numpy.float32(ObserverDistance), numpy.int32(BlackBody), numpy.int32(AngularCamera)) CLLaunch.wait() elapsed = time.time()-start_time print("Elapsed %f: " % elapsed) cl.enqueue_copy(queue,zImage,zImageCL).wait() cl.enqueue_copy(queue,fImage,fImageCL).wait() cl.enqueue_copy(queue,Trajectories,TrajectoriesCL).wait() cl.enqueue_copy(queue,IdLast,IdLastCL).wait() print(zImage.max()) print(fImage.max()) zImageCL.release() fImageCL.release() TrajectoriesCL.release() IdLastCL.release() return(elapsed) if __name__=='__main__': GpuStyle = 'OpenCL' Mass = 1. # Internal Radius 3 times de Schwarzschild Radius InternalRadius=6.*Mass # ExternalRadius=12. # ObserverDistance=100. # Angle with normal to disc 10 degrees Angle = numpy.pi/180.*(90.-10.) # Radiation of disc : BlackBody or Monochromatic BlackBody = True # Camera : Angular Camera or plate with the size of camera AngularCamera = True # Size of image Size=256 # Variable Type VariableType='FP32' # ? q=-2 # Method of resolution Method='Trajecto' HowToUse='%s -h [Help] -b [BlackBodyEmission] -c [AngularCamera] -s -m -i -x -o -a -d -g -t -v ' try: opts, args = getopt.getopt(sys.argv[1:],"hbcs:m:i:x:o:a:d:g:v:t:",["blackbody","camera","size=","mass=","internal=","external=","observer=","angle=","device=","gpustyle=","variabletype=","method="]) except getopt.GetoptError: print(HowToUse % sys.argv[0]) sys.exit(2) # List of Devices Devices=[] Alu={} for opt, arg in opts: if opt == '-h': print(HowToUse % sys.argv[0]) print("\nInformations about devices detected under OpenCL API:") # For PyOpenCL import try: import pyopencl as cl Id=0 for platform in cl.get_platforms(): for device in platform.get_devices(): #deviceType=cl.device_type.to_string(device.type) deviceType="xPU" print("Device #%i from %s of type %s : %s" % (Id,platform.vendor.lstrip(),deviceType,device.name.lstrip())) Id=Id+1 except: print("Your platform does not seem to support OpenCL") print("\nInformations about devices detected under CUDA API:") # For PyCUDA import try: import pycuda.driver as cuda cuda.init() for Id in range(cuda.Device.count()): device=cuda.Device(Id) print("Device #%i of type GPU : %s" % (Id,device.name())) print except: print("Your platform does not seem to support CUDA") sys.exit() elif opt in ("-d", "--device"): # Devices.append(int(arg)) Device=int(arg) elif opt in ("-g", "--gpustyle"): GpuStyle = arg elif opt in ("-v", "--variabletype"): VariableType = arg elif opt in ("-s", "--size"): Size = int(arg) elif opt in ("-m", "--mass"): Mass = float(arg) elif opt in ("-i", "--internal"): InternalRadius = float(arg) elif opt in ("-e", "--external"): ExternalRadius = float(arg) elif opt in ("-o", "--observer"): ObserverDistance = float(arg) elif opt in ("-a", "--angle"): Angle = numpy.pi/180.*(90.-float(arg)) elif opt in ("-b", "--blackbody"): BlackBody = True elif opt in ("-c", "--camera"): AngularCamera = True elif opt in ("-t", "--method"): Method = arg print("Device Identification selected : %s" % Device) print("GpuStyle used : %s" % GpuStyle) print("VariableType : %s" % VariableType) print("Size : %i" % Size) print("Mass : %f" % Mass) print("Internal Radius : %f" % InternalRadius) print("External Radius : %f" % ExternalRadius) print("Observer Distance : %f" % ObserverDistance) print("Angle with normal of (in radians) : %f" % Angle) print("Black Body Disc Emission (monochromatic instead) : %s" % BlackBody) print("Angular Camera (dimension of object instead) : %s" % AngularCamera) print("Method of resolution : %s" % Method) if GpuStyle=='CUDA': print("\nSelection of CUDA device") try: # For PyCUDA import import pycuda.driver as cuda cuda.init() for Id in range(cuda.Device.count()): device=cuda.Device(Id) print("Device #%i of type GPU : %s" % (Id,device.name())) if Id in Devices: Alu[Id]='GPU' except ImportError: print("Platform does not seem to support CUDA") if GpuStyle=='OpenCL': print("\nSelection of OpenCL device") try: # For PyOpenCL import import pyopencl as cl Id=0 for platform in cl.get_platforms(): for device in platform.get_devices(): #deviceType=cl.device_type.to_string(device.type) deviceType="xPU" print("Device #%i from %s of type %s : %s" % (Id,platform.vendor.lstrip().rstrip(),deviceType,device.name.lstrip().rstrip())) if Id in Devices: # Set the Alu as detected Device Type Alu[Id]=deviceType Id=Id+1 except ImportError: print("Platform does not seem to support OpenCL") # print(Devices,Alu) # MyImage=numpy.where(numpy.random.zeros(Size,Size)>0,1,-1).astype(numpy.float32) TrackPoints=2048 zImage=numpy.zeros((Size,Size),dtype=numpy.float32) fImage=numpy.zeros((Size,Size),dtype=numpy.float32) InputCL={} InputCL['Device']=Device InputCL['GpuStyle']=GpuStyle InputCL['VariableType']=VariableType InputCL['Size']=Size InputCL['Mass']=Mass InputCL['InternalRadius']=InternalRadius InputCL['ExternalRadius']=ExternalRadius InputCL['ObserverDistance']=ObserverDistance InputCL['Angle']=Angle InputCL['BlackBody']=BlackBody InputCL['AngularCamera']=AngularCamera InputCL['Method']=Method InputCL['TrackPoints']=TrackPoints duration=BlackHoleCL(zImage,fImage,InputCL) ImageOutput(zImage,"TrouNoirZ_%s" % Method) ImageOutput(fImage,"TrouNoirF_%s" % Method)