Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python

#

# TrouNoir model using PyOpenCL 

#

# CC BY-NC-SA 2019 : <emmanuel.quemener@ens-lyon.fr> 

#

# 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

import time,string

from numpy.random import randint as nprnd

import sys

import getopt

import matplotlib.pyplot as plt

from socket import gethostname

def DictionariesAPI():

    PhysicsList={'Einstein':0,'Newton':1}

    return(PhysicsList)

BlobOpenCL= """

#define PI (float)3.14159265359

#define nbr 256

#define EINSTEIN 0

#define NEWTON 1

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;

}

#if PHYSICS == NEWTON

float g(float u,float m,float b)

{

  return (-u);

}

#else

float g(float u,float m,float b)

{

  return (3.e0f*m/b*pow(u,2)-u);

}

#endif

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=exp(q*log(r/m))*pow(rf,4)*b*db*h;

  return(flx);

}

float spectre_cn(float rf32,float b32,float db32,

                 float h32,float r32,float m32,float bss32)

{

#define MYFLOAT double

  MYFLOAT rf=(MYFLOAT)(rf32);

  MYFLOAT b=(MYFLOAT)(b32);

  MYFLOAT db=(MYFLOAT)(db32);

  MYFLOAT h=(MYFLOAT)(h32);

  MYFLOAT r=(MYFLOAT)(r32);

  MYFLOAT m=(MYFLOAT)(m32);

  MYFLOAT bss=(MYFLOAT)(bss32);

  MYFLOAT flx;

  MYFLOAT nu_rec,nu_em,qu,temp_em,flux_int;

  int fi,posfreq;

#define planck 6.62e-34

#define k 1.38e-23

#define c2 9.e16

#define temp 3.e7

#define m_point 1.

#define lplanck (log(6.62)-34.*log(10.))

#define lk (log(1.38)-23.*log(10.))

#define lc2 (log(9.)+16.*log(10.))

  MYFLOAT v=1.-3./r;

  qu=1./sqrt((1.-3./r)*r)*(sqrt(r)-sqrt(6.)+sqrt(3.)/2.*log((sqrt(r)+sqrt(3.))/(sqrt(r)-sqrt(3.))* 0.17157287525380988 ));

  temp_em=temp*sqrt(m)*exp(0.25*log(m_point)-0.75*log(r)-0.125*log(v)+0.25*log(fabs(qu)));

  flux_int=0.;

  flx=0.;

  for (fi=0;fi<nbr;fi++)

    {

      nu_em=bss*(MYFLOAT)fi/(MYFLOAT)nbr;

      nu_rec=nu_em*rf;

      posfreq=(int)(nu_rec*(MYFLOAT)nbr/bss);

      if ((posfreq>0)&&(posfreq<nbr))

        {

          // Initial version

          // flux_int=2.*planck/c2*pow(nu_em,3)/(exp(planck*nu_em/(k*temp_em))-1.);

          // Version with log used

          //flux_int=2.*exp(lplanck-lc2+3.*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.);

          // flux_int*=pow(rf,3)*b*db*h;

          //flux_int*=exp(3.*log(rf))*b*db*h;

          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;

          flx+=flux_int;

        }

    }

  return((float)(flx));

}

void impact(float phi,float r,float b,float tho,float m,

            float *zp,float *fp,

            int q,float db,

            float h,int raie)

{

  float flx,rf,bss;

  rf=decalage_spectral(r,b,phi,tho,m);

  if (raie==0)

    {

      bss=1.e19;

      flx=spectre_cn(rf,b,db,h,r,m,bss);

    }

  else

    { 

      bss=2.;

      flx=spectre(rf,q,b,db,h,r,m,bss);

    }

  *zp=1./rf;

  *fp=flx;

}

__kernel void EachPixel(__global float *zImage,__global float *fImage,

                        float Mass,float InternalRadius,

                        float ExternalRadius,float Angle,

                        int Line)

{

   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, exit on hit

  float m,rs,ri,re,tho;

  int q,raie;

  m=Mass;

  rs=2.*m;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  q=-2;

  raie=Line;

  float d,bmx,db,b,h;

  float rp[2048];

  float phi,thi,phd,php,nr,r;

  int nh;

  float zp,fp;

  // Autosize for image

  bmx=1.25*re;

  b=0.;

  h=4.e0f*PI/(float)2048;

  // 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));

  float up,vp,pp,us,vs,ps;

  // impact parameter

  b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx;

  // step of impact parameter;

//  db=bmx/(float)(sizex/2);

  db=bmx/(float)(sizex);

  up=0.;

  vp=1.;

  pp=0.;

  nh=0;

  rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);

  rp[nh]=fabs(b/us);

  int ExitOnImpact=0;

  do

  {

     nh++;

     pp=ps;

     up=us;

     vp=vs;

     rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);          

     rp[nh]=fabs(b/us);

     ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rp[nh]>ri)&&(rp[nh]<re)?1:0;          

  } while ((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,

                    int Line)

{

   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;

  int q,raie;

  m=Mass;

  rs=2.*m;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  q=-2;

  raie=Line;

  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 external radius

  bmx=1.25*re;

  // Angular step of integration

  h=4.e0f*PI/(float)TrackPoints;

  // Step of Impact Parameter

  db=bmx/(2.e0*(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<ImpactParameter)

  {

    do

    {

      php=phd+(float)HalfLap*PI;

      nr=php/h;

      ni=(int)nr;

      if (ni<IdLast[bi])

      {

        r=(Trajectories[bi*TrackPoints+ni+1]-Trajectories[bi*TrackPoints+ni])*(nr-ni*1.)+Trajectories[bi*TrackPoints+ni];

      }

      else

      {

        r=Trajectories[bi*TrackPoints+ni];

      }

      if ((r<=re)&&(r>=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,

                     int Line)

{

   // 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;

  int q,raie;

  m=Mass;

  rs=2.*m;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  raie=Line;

  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/(2.e0*(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<IdLast[bi])

     {

        r=(Trajectories[bi*TrackPoints+ni+1]-Trajectories[bi*TrackPoints+ni])*(nr-ni*1.)+Trajectories[bi*TrackPoints+ni];

     }

     else

     {

        r=Trajectories[bi*TrackPoints+ni];

     }

     if ((r<=re)&&(r>=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,

                         int Line)

{

  // 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;

  int raie,q;

  m=Mass;

  rs=2.*m;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  q=-2;

  raie=Line;

  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;

  up=0.;

  vp=1.;

  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):

#     from PIL import Image

#     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 ImageOutput(sigma,prefix,Colors):

    import matplotlib.pyplot as plt

    if Colors == 'Red2Yellow':

        plt.imsave("%s.png" % prefix, sigma, cmap='afmhot')

    else:

        plt.imsave("%s.png" % prefix, sigma, cmap='Greys_r')

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']

    Method=InputCL['Method']

    TrackPoints=InputCL['TrackPoints']

    Physics=InputCL['Physics']

    PhysicsList=DictionariesAPI()

    if InputCL['BlackBody']:

        # Spectrum is Black Body one

        Line=0

    else:

        # Spectrum is Monochromatic Line one

        Line=1

    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()

    BuildOptions="-cl-mad-enable -DPHYSICS=%i " % (PhysicsList[Physics])

    BlackHoleCL = cl.Program(ctx,BlobOpenCL).build(options = BuildOptions)

    # Je recupere les flag possibles pour les buffers

    mf = cl.mem_flags

    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)

    IdLastCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=IdLast)

    start_time=time.time() 

    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.int32(Line))

        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.int32(Line))

        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.int32(Line))

        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.int32(Line))

        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.int32(Line))

        CLLaunch.wait()

    elapsed = time.time()-start_time

    print("\nElapsed Time : %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()

    zMaxPosition=numpy.where(zImage[:,:]==zImage.max())

    fMaxPosition=numpy.where(fImage[:,:]==fImage.max())

    print("Z max @(%i,%i) : %f" % (zMaxPosition[0],zMaxPosition[1],zImage.max()))

    print("Flux max @(%i,%i) : %f\n" % (fMaxPosition[0],fMaxPosition[1],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.

    #

    # Angle with normal to disc 10 degrees

    Angle = numpy.pi/180.*(90.-10.)

    # Radiation of disc : BlackBody or Monochromatic

    BlackBody = False

    # Size of image

    Size=256

    # Variable Type

    VariableType='FP32'

    # ?

    q=-2

    # Method of resolution

    Method='TrajectoPixel'

    # Colors for output image

    Colors='Greyscale'

    # Physics

    Physics='Einstein'

    # No output as image

    NoImage = False

    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>'

    try:

        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="])

    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 ("-a", "--angle"):

            Angle = numpy.pi/180.*(90.-float(arg))

        elif opt in ("-b", "--blackbody"):

            BlackBody = True

        elif opt in ("-n", "--noimage"):

            NoImage = True

        elif opt in ("-t", "--method"):

            Method = arg

        elif opt in ("-c", "--colors"):

            Colors = arg

        elif opt in ("-p", "--physics"):

            Physics = 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("Angle with normal of (in radians) : %f" % Angle)

    print("Black Body Disc Emission (monochromatic instead) : %s" % BlackBody)

    print("Method of resolution : %s" % Method)

    print("Colors for output images : %s" % Colors)

    print("Physics used for Trajectories : %s" % Physics)

    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['Angle']=Angle

    InputCL['BlackBody']=BlackBody

    InputCL['Method']=Method

    InputCL['TrackPoints']=TrackPoints

    InputCL['Physics']=Physics

    duration=BlackHoleCL(zImage,fImage,InputCL)

    Hostname=gethostname()

    Date=time.strftime("%Y%m%d_%H%M%S")

    ImageInfo="%s_Device%i_%s_%s" % (Method,Device,Hostname,Date)

    if not NoImage:

        ImageOutput(zImage,"TrouNoirZ_%s" % ImageInfo,Colors)

        ImageOutput(fImage,"TrouNoirF_%s" % ImageInfo,Colors)