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