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

# If crash on OpenCL Intel implementation, add following options and force

#export PYOPENCL_COMPILER_OUTPUT=1

#export CL_CONFIG_USE_VECTORIZER=True

#export CL_CONFIG_CPU_VECTORIZER_MODE=16

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)

#

# Blank space below to simplify debugging on OpenCL code

#

BlobOpenCL= """

#define PI (float)3.14159265359e0f

#define nbr 256

#define EINSTEIN 0

#define NEWTON 1

#ifdef SETTRACKPOINTS

#define TRACKPOINTS SETTRACKPOINTS

#else

#define TRACKPOINTS 2048

#endif

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

  c2=h*f(vp+c1/2.e0f);

  c3=h*f(vp+c2);

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

  d1=h*g(up+d0/2.e0f,m,b);

  d2=h*g(up+d1/2.e0f,m,b);

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

  *us=up+(c0+2.e0f*c1+2.e0f*c2+c3)/6.e0f;

  *vs=vp+(d0+2.e0f*d1+2.e0f*d2+d3)/6.e0f;

}

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.e0f*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-34f

#define k 1.38e-23f

#define c2 9.e16f

#define temp 3.e7f

#define m_point 1.e0f

#define lplanck (log(6.62e0f)-34.e0f*log(10.e0f))

#define lk (log(1.38e0f)-23.e0f*log(10.e0f))

#define lc2 (log(9.e0f)+16.e0f*log(10.e0f))

  MYFLOAT v=1.e0f-3.e0f/r;

  qu=1.e0f/sqrt((1.e0f-3.e0f/r)*r)*(sqrt(r)-sqrt(6.e0f)+sqrt(3.e0f)/2.e0f*log((sqrt(r)+sqrt(3.e0f))/(sqrt(r)-sqrt(3.e0f))* 0.17157287525380988e0f ));

  temp_em=temp*sqrt(m)*exp(0.25e0f*log(m_point)-0.75e0f*log(r)-0.125e0f*log(v)+0.25e0f*log(fabs(qu)));

  flux_int=0.e0f;

  flx=0.e0f;

  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.e0f*log(rf))*b*db*h;

          flux_int=2.e0f*exp(lplanck-lc2+3.e0f*log(nu_em))/(exp(exp(lplanck-lk+log(nu_em/temp_em)))-1.e0f)*exp(3.e0f*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.e19f;

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

    }

  else

    {

      bss=2.e0f;

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

    }

  *zp=1.e0f/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.e0f*m;

  ri=InternalRadius;

  re=ExternalRadius;

  tho=Angle;

  q=-2;

  raie=Line;

  float bmx,db,b,h;

  float rp0,rps;

  float phi,phd;

  uint nh=0;

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

  // Autosize for image

  bmx=1.25e0f*re;

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

  // set origin as center of image

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

  float y=(float)yi-(float)(sizey/2)+(float)5.e-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.e0f;

  vp=1.e0f;

  pp=0.e0f;

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

     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)&&(nh<TRACKPOINTS));

  if (ExitOnImpact==1) {

     impact(phi,rps,b,tho,m,&zp,&fp,q,db,h,raie);

  }

  else

  {

     zp=0.e0f;

     fp=0.e0f;

  }

  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,

                    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,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.e0f,fp=0.e0f;

  // 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)5.e-1f;

  float y=(float)yi-(float)(sizey/2)+(float)5.e-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));

  }

  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,

                     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,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.e0f,fp=0.e0f;

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

  barrier(CLK_GLOBAL_MEM_FENCE);

}

__kernel void Trajectory(__global float *Trajectories,__global int *IdLast,

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

  barrier(CLK_GLOBAL_MEM_FENCE);

}

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

                         float Mass,float InternalRadius,

                         float ExternalRadius,float Angle,

                         int Line)

{

   // Integer Impact Parameter ID

   uint bi=(uint)get_global_id(0);

   // Integer Impact Parameter Size (half of image)

   uint bmaxi=(uint)get_global_size(0);

  private float Trajectory[TRACKPOINTS];

  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;

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

  Trajectory[nh]=bvus;

  do

  {

     nh++;

     pp=ps;

     up=us;

     vp=vs;

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

     bvus=(float)fabs(b/us);

     Trajectory[nh]=bvus;

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

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

     Trajectory[i]=0.e0f;

  }

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

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

  {

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

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

     uint 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.e0f)+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;