Centre Blaise Pascal » Bench4GPU

#!/usr/bin/env python

#

#

#

# 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

# GPUfication in CUDA under Python in august 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

BlobOpenCL= """

#define nbr 256

#define EINSTEIN 0

#endif

void calcul(float *us,float *vs,float up,float vp,

{

  float c0,c1,c2,c3,d0,d1,d2,d3;

  c0=h*f(vp);

  c3=h*f(vp+c2);

  d0=h*g(up,m,b);

  d3=h*g(up+d2,m,b);

}

void rungekutta(float *ps,float *us,float *vs,

{

  calcul(us,vs,up,vp,h,m,b);

  *ps=pp+h;

}

float decalage_spectral(float r,float b,float phi,

{

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

//  flx=exp(q*log(r/m))*pow(rf,4)*b*db*h;

  return(flx);

}

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

{

#define MYFLOAT float

  MYFLOAT rf=(MYFLOAT)(rf32);

  MYFLOAT nu_rec,nu_em,qu,temp_em,flux_int;

  int fi,posfreq;

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

    {

      nu_rec=nu_em*rf;

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

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

          // Initial version

          // Version with log used

    }

  return((float)(flx));

}

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

{

  float flx,rf,bss;

  if (raie==0)

    {

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

    }

  else

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

    }

  *fp=flx;

}

  raie=Line;

  float bmx,db,b,h;

  float phi,phd;

  uint nh=0;

  float zp=0.e0f,fp=0.e0f;

  barrier(CLK_GLOBAL_MEM_FENCE);

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

}

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

  float bmx,db,b,h;

  float phi,phd,php,nr,r;

  // Autosize for image, 25% greater than external radius

  // Angular step of integration

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

  // Step of Impact Parameter

  // set origin as center of image

  // angle extracted from cylindric symmetry

  phi=atanp(x,y);

      if (ni<IdLast[bi])

      {

      }

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

                     float ExternalRadius,float Angle,

                     int Line)

{

   int bi=get_global_id(0);

   // Integer points on circle

   int i=get_global_id(1);

  float bmx,db,b,h;

  float phi,phd;

  // Autosize for image

  // Angular step of integration

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

  b=(float)bi/(float)bmaxi*bmx;

  db=bmx/(2.e0f*(float)bmaxi);

  phd=atanp(cos(phi)*sin(tho),cos(tho));

  int yi=(int)((float)bi*sin(phi))+bmaxi;

  int xi=(int)((float)bi*cos(phi))+bmaxi;

     if (ni<IdLast[bi])

     {

     }

     else

     {

     }

     if ((r<=re)&&(r>=ri))

     {

     }

     HalfLap++;

  } while ((HalfLap<=2)&&(ExitOnImpact==0));

  zImage[yi+2*bmaxi*xi]=zp;

  barrier(CLK_GLOBAL_MEM_FENCE);

                         float ExternalRadius,float Angle,

                         int Line)

{

  int bi=get_global_id(0);

  // Integer Impact Parameter Size (half of image)

  int bmaxi=get_global_size(0);

  float m,rs,re;

  m=Mass;

  re=ExternalRadius;

  float bmx,b,h;

  int nh;

  // Autosize for image

  // Angular step of integration

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

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

  nh=0;

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

  // b versus us

  float bvus=fabs(b/us);

  float bvus0=bvus;

  IdLast[bi]=nh;

  barrier(CLK_GLOBAL_MEM_FENCE);

}

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

                         float ExternalRadius,float Angle,

                         int Line)

{

  float m,rs,ri,re,tho;

  int raie,q;

  m=Mass;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  raie=Line;

  float bmx,db,b,h;

  // Autosize for image

  bmx=1.25e0f*re;

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

  nh=0;

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

  // b versus us

  float bvus=fabs(b/us);

  float bvus0=bvus;

     up=us;

     vp=vs;

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

     Trajectory[nh]=bvus;

  } while ((bvus>=rs)&&(bvus<=bvus0));

  for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {

     Trajectory[i]=0.e0f;

  }

  {

     float phd=atanp(cos(phi)*sin(tho),cos(tho));

     uint yi=(uint)((float)bi*sin(phi)+bmaxi);

     uint xi=(uint)((float)bi*cos(phi)+bmaxi);

     float php,nr,r;

     do

        if (ni<nh)

        {

        }

        else

        {

           r=Trajectory[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 Original(__global float *zImage,__global float *fImage,

   int raie,q;

   m=Mass;

   ri=InternalRadius;

   re=ExternalRadius;

   tho=Angle;

   raie=Line;

   float bmx,db,b,h;

   // Autosize for image

   bmx=1.25e0f*re;

   // Angular step of integration

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

   for (int bi=0;bi<bmaxi;bi++)

   {

      // impact parameter

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

      nh=0;

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

      // b versus us

      float bvus=fabs(b/us);

      float bvus0=bvus;

      for (int i=0;i<imx;i++)

      {

         float phi=2.e0f*PI/(float)imx*(float)i;

         float phd=atanp(cos(phi)*sin(tho),cos(tho));

         uint yi=(uint)((float)bi*sin(phi)+bmaxi);

         uint xi=(uint)((float)bi*cos(phi)+bmaxi);

         float php,nr,r;

         do

            if (ni<nh)

            {

            }

            else

            {

               r=Trajectory[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);

}

"""

#endif

__device__ void calcul(float *us,float *vs,float up,float vp,

{

  float c0,c1,c2,c3,d0,d1,d2,d3;

}

__device__ void rungekutta(float *ps,float *us,float *vs,

{

  calcul(us,vs,up,vp,h,m,b);

  *ps=pp+h;

}

__device__ float decalage_spectral(float r,float b,float phi,

{

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

}

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

{

#define MYFLOAT float

  MYFLOAT rf=(MYFLOAT)(rf32);

      nu_rec=nu_em*rf;

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

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

          // Initial version

          // Version with log used

    }

  return((float)(flx));

}

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

{

  float flx,rf,bss;

      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;

  up=0.;

  vp=1.;

  pp=0.;

  nh=0;

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

  rps=fabs(b/us);

  rp0=rps;

     up=us;

     vp=vs;

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

     rps=fabs(b/us);

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

  }

  else

  {

  }

  __syncthreads();

  int q,raie;

  m=Mass;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  int nh;

  float zp=0,fp=0;

  // Autosize for image, 25% greater than external radius

  // Angular step of integration

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

  // Step of Impact Parameter

  // set origin as center of image

  float x=(float)xi-(float)(sizex/2)+(float)5e-1f;

      if (ni<IdLast[bi])

      {

      }

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

                       float ExternalRadius,float Angle,

                       int Line)

{

   int bi=blockIdx.x*blockDim.x+threadIdx.x;

   // Integer points on circle

   int i=blockIdx.y*blockDim.y+threadIdx.y;

  int q,raie;

  m=Mass;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  float zp=0,fp=0;

  // Autosize for image

  // Angular step of integration

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

  // impact parameter

  b=(float)bi/(float)bmaxi*bmx;

  phd=atanp(cos(phi)*sin(tho),cos(tho));

  int yi=(int)((float)bi*sin(phi))+bmaxi;

  int xi=(int)((float)bi*cos(phi))+bmaxi;

     if (ni<IdLast[bi])

     {

     }

     else

     {

     }

     if ((r<=re)&&(r>=ri))

     {

     }

     HalfLap++;

  } while ((HalfLap<=2)&&(ExitOnImpact==0));

  zImage[yi+2*bmaxi*xi]=zp;

}

                           float ExternalRadius,float Angle,

                           int Line)

{

  int bi=blockIdx.x*blockDim.x+threadIdx.x;

  // Integer Impact Parameter Size (half of image)

  int bmaxi=gridDim.x*blockDim.x;

  int raie,q;

  m=Mass;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  float zp,fp;

  // Autosize for image

  // Angular step of integration

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

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

  nh=0;

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

  // b versus us

  float bvus=fabs(b/us);

  float bvus0=bvus;

                           float ExternalRadius,float Angle,

                           int Line)

{

  int bi=blockIdx.x*blockDim.x+threadIdx.x;

  // Integer Impact Parameter Size (half of image)

  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;

        {

           r=Trajectory[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));

   fImage[yi+2*bmaxi*xi]=fp;

  }

}

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

   int raie,q;

   m=Mass;

   ri=InternalRadius;

   re=ExternalRadius;

   tho=Angle;

   // Angular step of integration

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

   for (int bi=0;bi<bmaxi;bi++)

   {

      // impact parameter

      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;

            {

               r=Trajectory[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;

      }

}

"""

    return(BlobCUDA)

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

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

    save_time = time.time()-start_time

    print("Save Time : %f" % save_time)

def BlackHoleCL(zImage,fImage,InputCL):

    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

    if HasXPU==False:

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

    ctx = cl.Context([XPU])

    queue = cl.CommandQueue(ctx,

			    properties=cl.command_queue_properties.PROFILING_ENABLE)

    BuildOptions="-DPHYSICS=%i -DSETTRACKPOINTS=%i " % (PhysicsList[Physics],InputCL['TrackPoints'])

        print('No extra options for Intel and Clover!')

    else:

        BuildOptions = BuildOptions+" -cl-mad-enable"

    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)

    if Method=='EachPixel':

        CLLaunch=BlackHoleCL.EachPixel(queue,(zImage.shape[0],zImage.shape[1]),

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

    cuda.memcpy_htod(IdLastCU, IdLast)

    if Method=='EachPixel':

        EachPixelCU(zImageCU,fImageCU,