black-hole-accretion.py

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

  rps=fabs(b/us);
  rp0=rps;

  int ExitOnImpact=0;

  do
  {
     nh++;
     pp=ps;
     up=us;
     vp=vs;
     rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
     rpp=rps;
     rps=fabs(b/us);
     ExitOnImpact = ((fmod(pp,PI)<fmod(phd,PI))&&(fmod(ps,PI)>fmod(phd,PI)))&&(rps>ri)&&(rps<re)?1:0;

  } while ((rps>=rs)&&(rps<=rp0)&&(ExitOnImpact==0));

  if (ExitOnImpact==1) {
     impact(phi,rpp,b,tho,m,&zp,&fp,q,db,h,raie);
  }
  else
  {
     zp=0.e0f;
     fp=0.e0f;
  }

  __syncthreads();

  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,
                      float Mass,float InternalRadius,
                      float ExternalRadius,float Angle,
                      int Line)
{
   uint xi=(uint)(blockIdx.x*blockDim.x+threadIdx.x);
   uint yi=(uint)(blockIdx.y*blockDim.y+threadIdx.y);
   uint sizex=(uint)gridDim.x*blockDim.x;
   uint sizey=(uint)gridDim.y*blockDim.y;

  // Perform trajectory for each pixel

  float m,ri,re,tho;
  int q,raie;

  m=Mass;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  q=-2;
  raie=Line;

  float bmx,db,b,h;
  float phi,phd,php,nr,r;
  float zp=0,fp=0;
  // Autosize for image, 25% greater than external radius
  bmx=1.25e0f*re;

  // Angular step of integration
  h=4.e0f*PI/(float)TRACKPOINTS;

  // Step of Impact Parameter
  db=bmx/(2.e0f*(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.e0f)+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,
                       float Mass,float InternalRadius,
                       float ExternalRadius,float Angle,
                       int Line)
{
   // Integer Impact Parameter ID
   int bi=blockIdx.x*blockDim.x+threadIdx.x;
   // Integer points on circle
   int i=blockIdx.y*blockDim.y+threadIdx.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,ri,re,tho;
  int q,raie;

  m=Mass;
  ri=InternalRadius;
  re=ExternalRadius;
  tho=Angle;
  raie=Line;

  float bmx,db,b,h;
  float phi,phd;
  float zp=0,fp=0;

  // Autosize for image
  bmx=1.25e0f*re;

  // Angular step of integration
  h=4.e0f*PI/(float)TRACKPOINTS;

  // impact parameter
  b=(float)bi/(float)bmaxi*bmx;
  db=bmx/(2.e0f*(float)bmaxi);

  phi=2.e0f*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.e0f)+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,
                           float Mass,float InternalRadius,
                           float ExternalRadius,float Angle,
                           int Line)
{
  // Integer Impact Parameter ID
  int bi=blockIdx.x*blockDim.x+threadIdx.x;
  // Integer Impact Parameter Size (half of image)
  int bmaxi=gridDim.x*blockDim.x;

  // Perform trajectory for each pixel

  float m,rs,re;

  m=Mass;
  rs=2.e0f*m;
  re=ExternalRadius;

  float bmx,b,h;
  int nh;

  // Autosize for image
  bmx=1.25e0f*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.e0f;
  vp=1.e0f;
  pp=0.e0f;
  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;

}

__global__ void EachCircle(float *zImage,float *fImage,
                           float Mass,float InternalRadius,
                           float ExternalRadius,float Angle,
                           int Line)
{
  // Integer Impact Parameter ID
  int bi=blockIdx.x*blockDim.x+threadIdx.x;

  // Integer Impact Parameter Size (half of image)
  int bmaxi=gridDim.x*blockDim.x;

  float Trajectory[2048];

  // 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 bmx,db,b,h;
  int nh;

  // Autosize for image
  bmx=1.25e0f*re;

  // Angular step of integration
  h=4.e0f*PI/(float)TRACKPOINTS;

  // impact parameter
  b=(float)bi/(float)bmaxi*bmx;
  db=bmx/(2.e0f*(float)bmaxi);

  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;
  Trajectory[nh]=bvus;

  do
  {
     nh++;
     pp=ps;
     up=us;
     vp=vs;
     rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
     bvus=fabs(b/us);
     Trajectory[nh]=bvus;

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

  int imx=(int)(16*bi);

  for (int i=0;i<imx;i++)
  {
     float zp=0,fp=0;
     float phi=2.*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);

     int HalfLap=0,ExitOnImpact=0,ni;
     float php,nr,r;

     do
     {
        php=phd+(float)HalfLap*PI;
        nr=php/h;
        ni=(int)nr;

        if (ni<nh)
        {
           r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
        }
        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));

   __syncthreads();

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

  }

}

__global__ void Original(float *zImage,float *fImage,
                         uint Size,float Mass,float InternalRadius,
                         float ExternalRadius,float Angle,
                         int Line)
{
   // Integer Impact Parameter Size (half of image)
   uint bmaxi=(uint)Size;

   float Trajectory[TRACKPOINTS];

   // Perform trajectory for each pixel

   float m,rs,ri,re,tho;
   int raie,q;

   m=Mass;
   rs=2.e0f*m;
   ri=InternalRadius;
   re=ExternalRadius;
   tho=Angle;
   q=-2;
   raie=Line;

   float bmx,db,b,h;
   int nh;

   // Autosize for image
   bmx=1.25e0f*re;

   // Angular step of integration
   h=4.e0f*PI/(float)TRACKPOINTS;

   // Integer Impact Parameter ID
   for (int bi=0;bi<bmaxi;bi++)
   {
      // impact parameter
      b=(float)bi/(float)bmaxi*bmx;
      db=bmx/(2.e0f*(float)bmaxi);

      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;
      Trajectory[nh]=bvus;

      do
      {
         nh++;
         pp=ps;
         up=us;
         vp=vs;
         rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b);
         bvus=fabs(b/us);
         Trajectory[nh]=bvus;

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

      for (uint i=(uint)nh+1;i<TRACKPOINTS;i++) {
         Trajectory[i]=0.e0f;
      }

      int imx=(int)(16*bi);

      for (int i=0;i<imx;i++)
      {
         float zp=0,fp=0;
         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);

         int HalfLap=0,ExitOnImpact=0,ni;
         float php,nr,r;

         do
         {
            php=phd+(float)HalfLap*PI;
            nr=php/h;
            ni=(int)nr;

            if (ni<nh)
            {
               r=(Trajectory[ni+1]-Trajectory[ni])*(nr-ni*1.)+Trajectory[ni];
            }
            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;

      }

   }

}
"""
    return BlobCUDA


# def ImageOutput(sigma,prefix,Colors):
#     import matplotlib.pyplot as plt
#     start_time=time.time()
#     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 image as %s.png file" % prefix)
#     print("Save Time : %f" % save_time)


def ImageOutput(sigma, prefix, Colors):
    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 BlackHoleCL(zImage, fImage, InputCL):
    Device = InputCL["Device"]
    Mass = InputCL["Mass"]
    InternalRadius = InputCL["InternalRadius"]
    ExternalRadius = InputCL["ExternalRadius"]
    Angle = InputCL["Angle"]
    Method = InputCL["Method"]
    TrackPoints = InputCL["TrackPoints"]
    Physics = InputCL["Physics"]
    NoImage = InputCL["NoImage"]
    TrackSave = InputCL["TrackSave"]

    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)

    # 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:
                PF4XPU = platform.name
                XPU = device
                print("CPU/GPU selected: ", device.name.lstrip())
                HasXPU = True
            Id += 1

    if not HasXPU:
        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
    )

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

    print("My Platform is ", PF4XPU)

    if (
        "Intel" in PF4XPU
        or "Experimental" in PF4XPU
        or "Clover" in PF4XPU
        or "Portable" in PF4XPU
    ):
        print("No extra options for Intel and Clover!")
    else:
        BuildOptions = BuildOptions + " -cl-mad-enable"

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

    # Je recupere les flag possibles pour les buffers
    mf = cl.mem_flags

    if Method == "TrajectoPixel" or Method == "TrajectoCircle":
        TrajectoriesCL = cl.Buffer(
            ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=Trajectories
        )
        IdLastCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=IdLast)

    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)

    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 == "Original":
        CLLaunch = BlackHoleCL.Original(
            queue,
            (1,),
            None,
            zImageCL,
            fImageCL,
            numpy.uint32(zImage.shape[0] / 2),
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
        )
        CLLaunch.wait()
    elif Method == "EachCircle":
        CLLaunch = BlackHoleCL.EachCircle(
            queue,
            (int(zImage.shape[0] / 2),),
            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.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
        )

        CLLaunch = BlackHoleCL.Circle(
            queue,
            (Trajectories.shape[0], int(zImage.shape[0] * 4)),
            None,
            TrajectoriesCL,
            IdLastCL,
            zImageCL,
            fImageCL,
            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.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.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
        )
        CLLaunch.wait()

    compute = time.time() - start_time

    cl.enqueue_copy(queue, zImage, zImageCL).wait()
    cl.enqueue_copy(queue, fImage, fImageCL).wait()
    if Method == "TrajectoPixel" or Method == "TrajectoCircle":
        cl.enqueue_copy(queue, Trajectories, TrajectoriesCL).wait()
        cl.enqueue_copy(queue, IdLast, IdLastCL).wait()
    elapsed = time.time() - start_time
    print("\nCompute Time : %f" % compute)
    print("Elapsed Time : %f\n" % elapsed)

    zMaxPosition = numpy.where(zImage[:, :] == zImage.max())
    fMaxPosition = numpy.where(fImage[:, :] == fImage.max())
    print(
        "Z max @(%f,%f) : %f"
        % (
            (
                1.0 * zMaxPosition[1][0] / zImage.shape[1] - 0.5,
                1.0 * zMaxPosition[0][0] / zImage.shape[0] - 0.5,
                zImage.max(),
            )
        )
    )
    print(
        "Flux max @(%f,%f) : %f"
        % (
            (
                1.0 * fMaxPosition[1][0] / fImage.shape[1] - 0.5,
                1.0 * fMaxPosition[0][0] / fImage.shape[0] - 0.5,
                fImage.max(),
            )
        )
    )
    zImageCL.release()
    fImageCL.release()

    if Method == "TrajectoPixel" or Method == "TrajectoCircle":
        if not NoImage:
            AngleStep = 4 * numpy.pi / TrackPoints
            Angles = numpy.arange(0.0, 4 * numpy.pi, AngleStep)
            Angles.shape = (1, TrackPoints)

            if TrackSave:
                # numpy.savetxt("TrouNoirTrajectories_%s.csv" % ImageInfo,
                #               numpy.transpose(numpy.concatenate((Angles,Trajectories),axis=0)),
                #               delimiter=' ', fmt='%.2e')
                numpy.savetxt(
                    "TrouNoirTrajectories.csv",
                    numpy.transpose(numpy.concatenate((Angles, Trajectories),
                        axis=0)),
                    delimiter=" ",
                    fmt="%.2e",
                )

        TrajectoriesCL.release()
        IdLastCL.release()

    return elapsed


def BlackHoleCUDA(zImage, fImage, InputCL):
    Device = InputCL["Device"]
    Mass = InputCL["Mass"]
    InternalRadius = InputCL["InternalRadius"]
    ExternalRadius = InputCL["ExternalRadius"]
    Angle = InputCL["Angle"]
    Method = InputCL["Method"]
    TrackPoints = InputCL["TrackPoints"]
    Physics = InputCL["Physics"]
    Threads = InputCL["Threads"]

    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 -DSETTRACKPOINTS=%i"
                % (PhysicsList[Physics], TrackPoints),
            ],
        )
        print("Compilation seems to be OK")
    except Exception:
        print("Compilation seems to break")

    EachPixelCU = mod.get_function("EachPixel")
    OriginalCU = mod.get_function("Original")
    EachCircleCU = mod.get_function("EachCircle")
    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=(int(zImage.shape[0] / Threads), int(zImage.shape[1] / Threads)),
            block=(Threads, Threads, 1),
        )
    elif Method == "EachCircle":
        EachCircleCU(
            zImageCU,
            fImageCU,
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
            grid=(int(zImage.shape[0] / Threads / 2), 1),
            block=(Threads, 1, 1),
        )
    elif Method == "Original":
        OriginalCU(
            zImageCU,
            fImageCU,
            numpy.uint32(zImage.shape[0] / 2),
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
            grid=(1, 1),
            block=(1, 1, 1),
        )
    elif Method == "TrajectoCircle":
        TrajectoryCU(
            TrajectoriesCU,
            IdLastCU,
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
            grid=(int(Trajectories.shape[0] / Threads), 1),
            block=(Threads, 1, 1),
        )

        CircleCU(
            TrajectoriesCU,
            IdLastCU,
            zImageCU,
            fImageCU,
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
            grid=(
                int(Trajectories.shape[0] / Threads),
                int(zImage.shape[0] * 4 / Threads),
            ),
            block=(Threads, Threads, 1),
        )
    else:
        # Default method: TrajectoPixel
        TrajectoryCU(
            TrajectoriesCU,
            IdLastCU,
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
            grid=(int(Trajectories.shape[0] / Threads), 1),
            block=(Threads, 1, 1),
        )

        PixelCU(
            zImageCU,
            fImageCU,
            TrajectoriesCU,
            IdLastCU,
            numpy.uint32(Trajectories.shape[0]),
            numpy.float32(Mass),
            numpy.float32(InternalRadius),
            numpy.float32(ExternalRadius),
            numpy.float32(Angle),
            numpy.int32(Line),
            grid=(int(zImage.shape[0] / Threads), int(zImage.shape[1] / Threads), 1),
            block=(Threads, Threads, 1),
        )

    Context.synchronize()

    compute = time.time() - start_time

    cuda.memcpy_dtoh(zImage, zImageCU)
    cuda.memcpy_dtoh(fImage, fImageCU)
    if Method == "TrajectoPixel" or Method == "TrajectoCircle":
        cuda.memcpy_dtoh(Trajectories, TrajectoriesCU)
    elapsed = time.time() - start_time
    print("\nCompute Time : %f" % compute)
    print("Elapsed Time : %f\n" % elapsed)

    zMaxPosition = numpy.where(zImage[:, :] == zImage.max())
    fMaxPosition = numpy.where(fImage[:, :] == fImage.max())
    print(
        "Z max @(%f,%f) : %f"
        % (
            (
                1.0 * zMaxPosition[1][0] / zImage.shape[1] - 0.5,
                1.0 * zMaxPosition[0][0] / zImage.shape[0] - 0.5,
                zImage.max(),
            )
        )
    )
    print(
        "Flux max @(%f,%f) : %f"
        % (
            (
                1.0 * fMaxPosition[1][0] / fImage.shape[1] - 0.5,
                1.0 * fMaxPosition[0][0] / fImage.shape[0] - 0.5,
                fImage.max(),
            )
        )
    )

    Context.pop()

    Context.detach()

    if Method == "TrajectoPixel" or Method == "TrajectoCircle":
        if not NoImage:
            AngleStep = 4 * numpy.pi / TrackPoints
            Angles = numpy.arange(0.0, 4 * numpy.pi, AngleStep)
            Angles.shape = (1, TrackPoints)

            # numpy.savetxt("TrouNoirTrajectories_%s.csv" % ImageInfo,
            #               numpy.transpose(numpy.concatenate((Angles,Trajectories),axis=0)),
            #               delimiter=' ', fmt='%.2e')
            numpy.savetxt(
                "TrouNoirTrajectories.csv",
                numpy.transpose(numpy.concatenate((Angles, Trajectories), axis=0)),
                delimiter=" ",
                fmt="%.2e",
            )

    return elapsed


if __name__ == "__main__":

    # Default device: first one!
    Device = 0
    # Default implementation: OpenCL, most versatile!
    GpuStyle = "OpenCL"
    Mass = 1.0
    # Internal Radius 3 times de Schwarzschild Radius
    InternalRadius = 6.0 * Mass
    #
    ExternalRadius = 12.0
    #