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

# 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

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

#define TRACKPOINTS 2048

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;

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

  return(flx);

}

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

}

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

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

  private int q,raie;

  m=Mass;

  rs=2.*m;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  q=-2;

  raie=Line;

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

  private float rp0,rpp,rps;

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

  private int nh;

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

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

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

  flx=exp(q*log(r/m)+4.*log(rf))*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;