Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python3

#define TRACKPOINTS 2048

  private float m,rs,ri,re,tho;

  private int q,raie;

  private float d,bmx,db,b,h;

  private float rp[TRACKPOINTS];

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

  private int nh;

  private float zp,fp;

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

def KernelCodeCuda():

    BlobCUDA= """

#define PI (float)3.14159265359

#define nbr 256

#define EINSTEIN 0

#define NEWTON 1

#define TRACKPOINTS 2048

__device__ float nothing(float x)

{

  return(x);

}

__device__ float atanp(float x,float y)

{

  float angle;

  angle=atan2(y,x);

  if (angle<0.e0f)

    {

      angle+=(float)2.e0f*PI;

    }

  return(angle);

}

__device__ float f(float v)

{

  return(v);

}

#if PHYSICS == NEWTON

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

{

  return (-u);

}

#else

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

{

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

}

#endif

__device__ 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.;

}

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

}

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

}

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

}

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

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

{

#define MYFLOAT float

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

}

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

}

__global__ void EachPixel(float *zImage,float *fImage,

                          float Mass,float InternalRadius,

                          float ExternalRadius,float Angle,

                          int Line)

{

   uint xi=(uint)blockIdx.x;

   uint yi=(uint)blockIdx.y;

   uint sizex=(uint)gridDim.x*blockDim.x;

   uint sizey=(uint)gridDim.y*blockDim.y;

   // 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[TRACKPOINTS];

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

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

  }

  zImage[yi+sizex*xi]=(float)zp;

  fImage[yi+sizex*xi]=(float)fp;

}

__global__ void Pixel(float *zImage,float *fImage,

                      float *Trajectories,int *IdLast,

                      uint ImpactParameter,uint TrackPoints,

                      float Mass,float InternalRadius,

                      float ExternalRadius,float Angle,

                      int Line)

{

  uint xi=(uint)blockIdx.x;

  uint yi=(uint)blockIdx.y;

  uint sizex=(uint)gridDim.x*blockDim.x;

  uint sizey=(uint)gridDim.y*blockDim.y;

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

  }

  zImage[yi+sizex*xi]=zp;

  fImage[yi+sizex*xi]=fp;

}

__global__ void Circle(float *Trajectories,int *IdLast,

                       float *zImage,float *fImage,

                       int TrackPoints,

                       float Mass,float InternalRadius,

                       float ExternalRadius,float Angle,

                       int Line)

{

   // Integer Impact Parameter ID 

   int bi=blockIdx.x;

   // Integer points on circle

   int i=blockIdx.y;

   // Integer Impact Parameter Size (half of image)

   int bmaxi=gridDim.x*blockDim.x;

   // Integer Points on circle

   int imx=gridDim.y*blockDim.y;

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

}

__global__ void Trajectory(float *Trajectories,

                           int *IdLast,int TrackPoints,

                           float Mass,float InternalRadius,

                           float ExternalRadius,float Angle,

                           int Line)

{

  // Integer Impact Parameter ID 

  int bi=blockIdx.x;

  // Integer Impact Parameter Size (half of image)

  int bmaxi=gridDim.x*blockDim.x;

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

}

"""

    return(BlobCUDA)

    Trajectories=numpy.zeros((int(InputCL['Size']/2),InputCL['TrackPoints']),dtype=numpy.float32)

    IdLast=numpy.zeros(int(InputCL['Size']/2),dtype=numpy.int32)

                print("CPU/GPU selected: ",device.name.lstrip())

        print("No XPU #%i found in all of %i devices, sorry..." % (Device,Id-1))

                                          zImage.shape[1]),None,

def BlackHoleCUDA(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((int(InputCL['Size']/2),InputCL['TrackPoints']),dtype=numpy.float32)

    IdLast=numpy.zeros(int(InputCL['Size']/2),dtype=numpy.int32)

    try:

        # For PyCUDA import

        import pycuda.driver as cuda

        from pycuda.compiler import SourceModule

        cuda.init()

        for Id in range(cuda.Device.count()):

            if Id==Device:

                XPU=cuda.Device(Id)

                print("GPU selected %s" % XPU.name())

        print

    except ImportError:

        print("Platform does not seem to support CUDA")

    Context=XPU.make_context()

    try:

        mod = SourceModule(KernelCodeCuda(),options=['--compiler-options','-DPHYSICS=%i' % (PhysicsList[Physics])])

        print("Compilation seems to be OK")

    except:

        print("Compilation seems to break")

    EachPixelCU=mod.get_function("EachPixel")

    TrajectoryCU=mod.get_function("Trajectory")

    PixelCU=mod.get_function("Pixel")

    CircleCU=mod.get_function("Circle")

    TrajectoriesCU = cuda.mem_alloc(Trajectories.size*Trajectories.dtype.itemsize)

    cuda.memcpy_htod(TrajectoriesCU, Trajectories)

    zImageCU = cuda.mem_alloc(zImage.size*zImage.dtype.itemsize)

    cuda.memcpy_htod(zImageCU, zImage)

    fImageCU = cuda.mem_alloc(fImage.size*fImage.dtype.itemsize)

    cuda.memcpy_htod(zImageCU, fImage)

    IdLastCU = cuda.mem_alloc(IdLast.size*IdLast.dtype.itemsize)

    cuda.memcpy_htod(IdLastCU, IdLast)

    start_time=time.time() 

    if Method=='EachPixel':

        EachPixelCU(zImageCU,fImageCU,

                    numpy.float32(Mass),

                    numpy.float32(InternalRadius),

                    numpy.float32(ExternalRadius),

                    numpy.float32(Angle),

                    numpy.int32(Line),

                    grid=(zImage.shape[0],zImage.shape[1]),block=(1,1,1))

    elif Method=='TrajectoCircle':

        TrajectoryCU(TrajectoriesCU,IdLastCU,

                     numpy.uint32(Trajectories.shape[1]),

                     numpy.float32(Mass),

                     numpy.float32(InternalRadius),

                     numpy.float32(ExternalRadius),

                     numpy.float32(Angle),

                     numpy.int32(Line),

                     grid=(Trajectories.shape[0],1),block=(1,1,1))

        CircleCU(TrajectoriesCU,IdLastCU,zImageCU,fImageCU,

                 numpy.uint32(Trajectories.shape[1]),

                 numpy.float32(Mass),

                 numpy.float32(InternalRadius),

                 numpy.float32(ExternalRadius),

                 numpy.float32(Angle),

                 numpy.int32(Line),

                 grid=(Trajectories.shape[0],zImage.shape[0]*4),block=(1,1,1))

    else:

        TrajectoryCU(TrajectoriesCU,IdLastCU,

                     numpy.uint32(Trajectories.shape[1]),

                     numpy.float32(Mass),

                     numpy.float32(InternalRadius),

                     numpy.float32(ExternalRadius),

                     numpy.float32(Angle),

                     numpy.int32(Line),

                     grid=(Trajectories.shape[0],1),block=(1,1,1))

        PixelCU(zImageCU,fImageCU,TrajectoriesCU,IdLastCU,

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

                grid=(zImage.shape[0],zImage.shape[1],1),block=(1,1,1))

    elapsed = time.time()-start_time

    print("\nElapsed Time : %f" % elapsed)

    cuda.memcpy_dtoh(zImage,zImageCU)

    cuda.memcpy_dtoh(fImage,fImageCU)

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

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

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

    print("Flux max @(%i,%i) : %f\n" % (fMaxPosition[0][0],fMaxPosition[1][0],fImage.max()))

    Context.detach()

    return(elapsed)

    if GpuStyle=='OpenCL':

        duration=BlackHoleCL(zImage,fImage,InputCL)

    else:

        duration=BlackHoleCUDA(zImage,fImage,InputCL)