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

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;fi<nbr;fi++)

    {

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

      nu_rec=nu_em*rf;

      posfreq=(int)(1./bss*nu_rec*(float)nbr);

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

        {

          //flux_int=2*planck/9.e16f*pow(nu_em,3)/(exp(planck*nu_em/k/temp_em)-1.)*pow(rf,3)*b*db*h;

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

          flux_int=2.*psc2*pow(nu_em,3)/(exp(psk*nu_em/temp_em)-1.)*pow(rf,3)*b;

          flx+=flux_int;

        }

    }

  //return(pow(rf,3)*b*temp_em);

  return((float)flx);

  return(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=2.;

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

    }

  else

    {

      bss=1.e6;

      flx=spectre_cn(rf,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,

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

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

  if (dist==1)

  {

     up=0.;

     vp=1.;

  }

  else

  {

     up=sin(thi);

     vp=cos(thi);

  }

  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,

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

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

                         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 <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -o <ObservatorDistance> -a <AngleAboveDisc> -d <DeviceId> -g <CUDA/OpenCL> -t <EachPixel/TrajectoCircle/TrajectoPixel> -v <FP32/FP64>'

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