Révision 203
| TrouNoir/TrouNoir.py (revision 203) | ||
|---|---|---|
| 2 | 2 |
# |
| 3 | 3 |
# TrouNoir model using PyOpenCL |
| 4 | 4 |
# |
| 5 |
# CC BY-NC-SA 2011 : <emmanuel.quemener@ens-lyon.fr>
|
|
| 5 |
# CC BY-NC-SA 2019 : <emmanuel.quemener@ens-lyon.fr>
|
|
| 6 | 6 |
# |
| 7 | 7 |
# Thanks to Andreas Klockner for PyOpenCL: |
| 8 | 8 |
# http://mathema.tician.de/software/pyopencl |
| 9 | 9 |
# |
| 10 |
# Interesting links: |
|
| 11 |
# http://viennacl.sourceforge.net/viennacl-documentation.html |
|
| 12 |
# http://enja.org/2011/02/22/adventures-in-pyopencl-part-1-getting-started-with-python/ |
|
| 10 |
# Original code programmed in Fortran 77 in mars 1994 |
|
| 11 |
# for Practical Work of Numerical Simulation |
|
| 12 |
# DEA (old Master2) in astrophysics and spatial techniques in Meudon |
|
| 13 |
# by Herve Aussel & Emmanuel Quemener |
|
| 14 |
# |
|
| 15 |
# Conversion in C done by Emmanuel Quemener in august 1997 |
|
| 16 |
# GPUfication in OpenCL under Python in july 2019 |
|
| 17 |
# |
|
| 18 |
# Thanks to : |
|
| 19 |
# |
|
| 20 |
# - Herve Aussel for his part of code of black body spectrum |
|
| 21 |
# - Didier Pelat for his help to perform this work |
|
| 22 |
# - Jean-Pierre Luminet for his article published in 1979 |
|
| 23 |
# - Numerical Recipies for Runge Kutta recipies |
|
| 24 |
# - Luc Blanchet for his disponibility about my questions in General Relativity |
|
| 25 |
# - Pierre Lena for his passion about science and vulgarisation |
|
| 13 | 26 |
|
| 14 | 27 |
import pyopencl as cl |
| 15 | 28 |
import numpy |
| ... | ... | |
| 20 | 33 |
import getopt |
| 21 | 34 |
import matplotlib.pyplot as plt |
| 22 | 35 |
|
| 36 |
|
|
| 37 |
|
|
| 38 |
|
|
| 39 |
|
|
| 40 |
|
|
| 41 |
|
|
| 42 |
|
|
| 23 | 43 |
KERNEL_CODE=string.Template("""
|
| 24 | 44 |
|
| 25 | 45 |
#define PI (float)3.14159265359 |
| ... | ... | |
| 31 | 51 |
|
| 32 | 52 |
angle=atan2(y,x); |
| 33 | 53 |
|
| 34 |
if (angle<0) |
|
| 54 |
if (angle<0.e0f)
|
|
| 35 | 55 |
{
|
| 36 | 56 |
angle+=(float)2.e0f*PI; |
| 37 | 57 |
} |
| ... | ... | |
| 46 | 66 |
|
| 47 | 67 |
float g(float u,float m,float b) |
| 48 | 68 |
{
|
| 49 |
return (3.*m/b*pow(u,2)-u); |
|
| 69 |
return (3.e0f*m/b*pow(u,2)-u);
|
|
| 50 | 70 |
} |
| 51 | 71 |
|
| 52 | 72 |
|
| ... | ... | |
| 91 | 111 |
return(flx); |
| 92 | 112 |
} |
| 93 | 113 |
|
| 94 |
float spectre_cn(float rf,float b,float db, |
|
| 95 |
float h,float r,float m,float bss) |
|
| 114 |
float spectre_cn(float rf,float b,float db,float h,float r,float m,float bss) |
|
| 96 | 115 |
{
|
| 97 | 116 |
double flx; |
| 98 | 117 |
double nu_rec,nu_em,qu,v,temp,temp_em,flux_int,m_point,planck,c,k,psc2,psk; |
| ... | ... | |
| 134 | 153 |
void impact(float phi,float r,float b,float tho,float m, |
| 135 | 154 |
float *zp,float *fp, |
| 136 | 155 |
int q,float db, |
| 137 |
float h,float bss,int raie)
|
|
| 156 |
float h,int raie) |
|
| 138 | 157 |
{
|
| 139 |
float flx,rf; |
|
| 158 |
float flx,rf,bss;
|
|
| 140 | 159 |
|
| 141 | 160 |
rf=decalage_spectral(r,b,phi,tho,m); |
| 142 | 161 |
|
| 143 | 162 |
if (raie==0) |
| 144 | 163 |
{
|
| 164 |
bss=2.; |
|
| 145 | 165 |
flx=spectre(rf,q,b,db,h,r,m,bss); |
| 146 | 166 |
} |
| 147 | 167 |
else |
| 148 | 168 |
{
|
| 169 |
bss=1.e6; |
|
| 149 | 170 |
flx=spectre_cn(rf,b,db,h,r,m,bss); |
| 150 | 171 |
} |
| 151 | 172 |
|
| ... | ... | |
| 154 | 175 |
|
| 155 | 176 |
} |
| 156 | 177 |
|
| 157 |
__kernel void EachPixel(__global float *zImage,__global float *fImage) |
|
| 178 |
__kernel void EachPixel(__global float *zImage,__global float *fImage, |
|
| 179 |
float Mass,float InternalRadius, |
|
| 180 |
float ExternalRadius,float Angle, |
|
| 181 |
float ObserverDistance, |
|
| 182 |
int BlackBody,int AngularCamera) |
|
| 158 | 183 |
{
|
| 159 | 184 |
uint xi=(uint)get_global_id(0); |
| 160 | 185 |
uint yi=(uint)get_global_id(1); |
| 161 | 186 |
uint sizex=(uint)get_global_size(0); |
| 162 | 187 |
uint sizey=(uint)get_global_size(1); |
| 163 | 188 |
|
| 164 |
// Perform trajectory for each pixel |
|
| 189 |
// Perform trajectory for each pixel, exit on hit
|
|
| 165 | 190 |
|
| 166 | 191 |
float m,rs,ri,re,tho,ro; |
| 167 |
int q; |
|
| 192 |
int q,dist,raie;
|
|
| 168 | 193 |
|
| 169 |
float bss,stp; |
|
| 170 |
int nmx,dim; |
|
| 194 |
m=Mass; |
|
| 195 |
rs=2.*m; |
|
| 196 |
ri=InternalRadius; |
|
| 197 |
re=ExternalRadius; |
|
| 198 |
ro=ObserverDistance; |
|
| 199 |
tho=Angle; |
|
| 200 |
q=-2; |
|
| 201 |
dist=AngularCamera; |
|
| 202 |
raie=BlackBody; |
|
| 203 |
|
|
| 171 | 204 |
float d,bmx,db,b,h; |
| 172 |
float up,vp,pp; |
|
| 173 |
float us,vs,ps; |
|
| 174 |
float rp[2000]; |
|
| 175 |
float zmx,fmx,zen,fen; |
|
| 176 |
float flux_tot,impc_tot; |
|
| 177 |
float phi,thi,thx,phd,php,nr,r; |
|
| 178 |
int ni,ii,i,imx,j,n,tst,dist,raie,pc,fcl,zcl; |
|
| 205 |
float rp[2048]; |
|
| 206 |
float phi,thi,phd,php,nr,r; |
|
| 179 | 207 |
int nh; |
| 180 |
int ExitOnImpact; |
|
| 181 | 208 |
float zp,fp; |
| 182 | 209 |
|
| 183 |
m=1.; |
|
| 184 |
rs=2.*m; |
|
| 185 |
ri=3.*rs; |
|
| 186 |
re=12.; |
|
| 187 |
ro=100.; |
|
| 188 |
//tho=PI/180.*80; |
|
| 189 |
tho=PI/180.*80.; |
|
| 190 |
q=-2; |
|
| 191 |
dist=0; |
|
| 192 |
raie=0; |
|
| 193 |
|
|
| 194 |
nmx=(float)sizex/2.; |
|
| 195 |
stp=0.5; |
|
| 196 | 210 |
// Autosize for image |
| 197 | 211 |
bmx=1.25*re; |
| 198 | 212 |
b=0.; |
| 199 |
thx=asin(bmx/ro); |
|
| 200 |
pc=0; |
|
| 201 | 213 |
|
| 202 |
if (raie==0) |
|
| 203 |
{
|
|
| 204 |
bss=2.; |
|
| 205 |
} |
|
| 206 |
else |
|
| 207 |
{
|
|
| 208 |
bss=1.e6; |
|
| 209 |
} |
|
| 214 |
h=4.e0f*PI/(float)2048; |
|
| 210 | 215 |
|
| 211 |
h=PI/5.e2f; |
|
| 212 |
|
|
| 213 | 216 |
// set origin as center of image |
| 214 | 217 |
float x=(float)xi-(float)(sizex/2)+(float)5e-1f; |
| 215 | 218 |
float y=(float)yi-(float)(sizey/2)+(float)5e-1f; |
| ... | ... | |
| 217 | 220 |
phi=atanp(x,y); |
| 218 | 221 |
phd=atanp(cos(phi)*sin(tho),cos(tho)); |
| 219 | 222 |
|
| 223 |
float up,vp,pp,us,vs,ps; |
|
| 224 |
|
|
| 225 |
// impact parameter |
|
| 226 |
b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx; |
|
| 227 |
// step of impact parameter; |
|
| 228 |
db=bmx/(float)(sizex/2); |
|
| 229 |
|
|
| 220 | 230 |
if (dist==1) |
| 221 | 231 |
{
|
| 222 |
db=bmx/nmx; |
|
| 223 |
// impact parameter |
|
| 224 |
b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx; |
|
| 225 | 232 |
up=0.; |
| 226 | 233 |
vp=1.; |
| 227 | 234 |
} |
| 228 | 235 |
else |
| 229 | 236 |
{
|
| 230 |
b=sqrt(x*x+y*y)*(float)2.e0f/(float)sizex*bmx; |
|
| 231 |
thi=asin(b/ro); |
|
| 232 |
db=bmx/nmx; |
|
| 233 | 237 |
up=sin(thi); |
| 234 | 238 |
vp=cos(thi); |
| 235 | 239 |
} |
| ... | ... | |
| 241 | 245 |
|
| 242 | 246 |
rp[nh]=fabs(b/us); |
| 243 | 247 |
|
| 244 |
ExitOnImpact=0; |
|
| 248 |
int ExitOnImpact=0;
|
|
| 245 | 249 |
|
| 246 | 250 |
do |
| 247 | 251 |
{
|
| ... | ... | |
| 256 | 260 |
} while ((rp[nh]>=rs)&&(rp[nh]<=rp[0])&&(ExitOnImpact==0)); |
| 257 | 261 |
|
| 258 | 262 |
if (ExitOnImpact==1) {
|
| 259 |
impact(phi,rp[nh-1],b,tho,m,&zp,&fp,q,db,h,bss,raie);
|
|
| 263 |
impact(phi,rp[nh-1],b,tho,m,&zp,&fp,q,db,h,raie); |
|
| 260 | 264 |
} |
| 261 | 265 |
else |
| 262 | 266 |
{
|
| ... | ... | |
| 264 | 268 |
fp=0.; |
| 265 | 269 |
} |
| 266 | 270 |
|
| 267 |
// |
|
| 271 |
barrier(CLK_GLOBAL_MEM_FENCE); |
|
| 272 |
|
|
| 268 | 273 |
zImage[yi+sizex*xi]=(float)zp; |
| 269 |
fImage[yi+sizex*xi]=(float)fp; |
|
| 274 |
fImage[yi+sizex*xi]=(float)fp; |
|
| 275 |
} |
|
| 270 | 276 |
|
| 277 |
__kernel void Pixel(__global float *zImage,__global float *fImage, |
|
| 278 |
__global float *Trajectories,__global int *IdLast, |
|
| 279 |
uint ImpactParameter,uint TrackPoints, |
|
| 280 |
float Mass,float InternalRadius, |
|
| 281 |
float ExternalRadius,float Angle, |
|
| 282 |
float ObserverDistance, |
|
| 283 |
int BlackBody,int AngularCamera) |
|
| 284 |
{
|
|
| 285 |
uint xi=(uint)get_global_id(0); |
|
| 286 |
uint yi=(uint)get_global_id(1); |
|
| 287 |
uint sizex=(uint)get_global_size(0); |
|
| 288 |
uint sizey=(uint)get_global_size(1); |
|
| 289 |
|
|
| 290 |
// Perform trajectory for each pixel |
|
| 291 |
|
|
| 292 |
float m,rs,ri,re,tho,ro; |
|
| 293 |
int q,dist,raie; |
|
| 294 |
|
|
| 295 |
m=Mass; |
|
| 296 |
rs=2.*m; |
|
| 297 |
ri=InternalRadius; |
|
| 298 |
re=ExternalRadius; |
|
| 299 |
ro=ObserverDistance; |
|
| 300 |
tho=Angle; |
|
| 301 |
q=-2; |
|
| 302 |
dist=AngularCamera; |
|
| 303 |
raie=BlackBody; |
|
| 304 |
|
|
| 305 |
float d,bmx,db,b,h; |
|
| 306 |
float phi,thi,phd,php,nr,r; |
|
| 307 |
int nh; |
|
| 308 |
float zp=0,fp=0; |
|
| 309 |
|
|
| 310 |
// Autosize for image, 25% greater than externa radius |
|
| 311 |
bmx=1.25*re; |
|
| 312 |
|
|
| 313 |
// Angular step of integration |
|
| 314 |
h=4.e0f*PI/(float)TrackPoints; |
|
| 315 |
|
|
| 316 |
// Step of Impact Parameter |
|
| 317 |
db=bmx/(float)ImpactParameter; |
|
| 318 |
|
|
| 319 |
// set origin as center of image |
|
| 320 |
float x=(float)xi-(float)(sizex/2)+(float)5e-1f; |
|
| 321 |
float y=(float)yi-(float)(sizey/2)+(float)5e-1f; |
|
| 322 |
|
|
| 323 |
// angle extracted from cylindric symmetry |
|
| 324 |
phi=atanp(x,y); |
|
| 325 |
phd=atanp(cos(phi)*sin(tho),cos(tho)); |
|
| 326 |
|
|
| 327 |
// Real Impact Parameter |
|
| 328 |
b=sqrt(x*x+y*y)*bmx/(float)ImpactParameter; |
|
| 329 |
|
|
| 330 |
// Integer Impact Parameter |
|
| 331 |
uint bi=(uint)sqrt(x*x+y*y); |
|
| 332 |
|
|
| 333 |
int HalfLap=0,ExitOnImpact=0,ni; |
|
| 334 |
|
|
| 335 |
if (bi<ImpactParameter) |
|
| 336 |
{
|
|
| 337 |
do |
|
| 338 |
{
|
|
| 339 |
php=phd+(float)HalfLap*PI; |
|
| 340 |
nr=php/h; |
|
| 341 |
ni=(int)nr; |
|
| 342 |
|
|
| 343 |
if (ni<IdLast[bi]) |
|
| 344 |
{
|
|
| 345 |
r=(Trajectories[bi*TrackPoints+ni+1]-Trajectories[bi*TrackPoints+ni])*(nr-ni*1.)+Trajectories[bi*TrackPoints+ni]; |
|
| 346 |
} |
|
| 347 |
else |
|
| 348 |
{
|
|
| 349 |
r=Trajectories[bi*TrackPoints+ni]; |
|
| 350 |
} |
|
| 351 |
|
|
| 352 |
if ((r<=re)&&(r>=ri)) |
|
| 353 |
{
|
|
| 354 |
ExitOnImpact=1; |
|
| 355 |
impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie); |
|
| 356 |
} |
|
| 357 |
|
|
| 358 |
HalfLap++; |
|
| 359 |
} while ((HalfLap<=2)&&(ExitOnImpact==0)); |
|
| 360 |
|
|
| 361 |
} |
|
| 362 |
|
|
| 271 | 363 |
barrier(CLK_GLOBAL_MEM_FENCE); |
| 364 |
|
|
| 365 |
zImage[yi+sizex*xi]=zp; |
|
| 366 |
fImage[yi+sizex*xi]=fp; |
|
| 367 |
} |
|
| 368 |
|
|
| 369 |
__kernel void Circle(__global float *Trajectories,__global int *IdLast, |
|
| 370 |
__global float *zImage,__global float *fImage, |
|
| 371 |
int TrackPoints, |
|
| 372 |
float Mass,float InternalRadius, |
|
| 373 |
float ExternalRadius,float Angle, |
|
| 374 |
float ObserverDistance, |
|
| 375 |
int BlackBody,int AngularCamera) |
|
| 376 |
{
|
|
| 377 |
// Integer Impact Parameter ID |
|
| 378 |
int bi=get_global_id(0); |
|
| 379 |
// Integer points on circle |
|
| 380 |
int i=get_global_id(1); |
|
| 381 |
// Integer Impact Parameter Size (half of image) |
|
| 382 |
int bmaxi=get_global_size(0); |
|
| 383 |
// Integer Points on circle |
|
| 384 |
int imx=get_global_size(1); |
|
| 385 |
|
|
| 386 |
// Perform trajectory for each pixel |
|
| 387 |
|
|
| 388 |
float m,rs,ri,re,tho,ro; |
|
| 389 |
int q,dist,raie; |
|
| 390 |
|
|
| 391 |
m=Mass; |
|
| 392 |
rs=2.*m; |
|
| 393 |
ri=InternalRadius; |
|
| 394 |
re=ExternalRadius; |
|
| 395 |
ro=ObserverDistance; |
|
| 396 |
tho=Angle; |
|
| 397 |
q=-2; |
|
| 398 |
dist=AngularCamera; |
|
| 399 |
raie=BlackBody; |
|
| 400 |
|
|
| 401 |
float bmx,db,b,h; |
|
| 402 |
float phi,thi,phd; |
|
| 403 |
int nh; |
|
| 404 |
float zp=0,fp=0; |
|
| 405 |
|
|
| 406 |
// Autosize for image |
|
| 407 |
bmx=1.25*re; |
|
| 408 |
|
|
| 409 |
// Angular step of integration |
|
| 410 |
h=4.e0f*PI/(float)TrackPoints; |
|
| 411 |
|
|
| 412 |
// impact parameter |
|
| 413 |
b=(float)bi/(float)bmaxi*bmx; |
|
| 414 |
db=bmx/(float)bmaxi; |
|
| 415 |
|
|
| 416 |
phi=2.*PI/(float)imx*(float)i; |
|
| 417 |
phd=atanp(cos(phi)*sin(tho),cos(tho)); |
|
| 418 |
int yi=(int)((float)bi*sin(phi))+bmaxi; |
|
| 419 |
int xi=(int)((float)bi*cos(phi))+bmaxi; |
|
| 420 |
|
|
| 421 |
int HalfLap=0,ExitOnImpact=0,ni; |
|
| 422 |
float php,nr,r; |
|
| 423 |
|
|
| 424 |
do |
|
| 425 |
{
|
|
| 426 |
php=phd+(float)HalfLap*PI; |
|
| 427 |
nr=php/h; |
|
| 428 |
ni=(int)nr; |
|
| 429 |
|
|
| 430 |
if (ni<IdLast[bi]) |
|
| 431 |
{
|
|
| 432 |
r=(Trajectories[bi*TrackPoints+ni+1]-Trajectories[bi*TrackPoints+ni])*(nr-ni*1.)+Trajectories[bi*TrackPoints+ni]; |
|
| 433 |
} |
|
| 434 |
else |
|
| 435 |
{
|
|
| 436 |
r=Trajectories[bi*TrackPoints+ni]; |
|
| 437 |
} |
|
| 438 |
|
|
| 439 |
if ((r<=re)&&(r>=ri)) |
|
| 440 |
{
|
|
| 441 |
ExitOnImpact=1; |
|
| 442 |
impact(phi,r,b,tho,m,&zp,&fp,q,db,h,raie); |
|
| 443 |
} |
|
| 444 |
|
|
| 445 |
HalfLap++; |
|
| 446 |
} while ((HalfLap<=2)&&(ExitOnImpact==0)); |
|
| 447 |
|
|
| 448 |
zImage[yi+2*bmaxi*xi]=zp; |
|
| 449 |
fImage[yi+2*bmaxi*xi]=fp; |
|
| 450 |
|
|
| 451 |
barrier(CLK_GLOBAL_MEM_FENCE); |
|
| 452 |
|
|
| 453 |
} |
|
| 454 |
|
|
| 455 |
__kernel void Trajectory(__global float *Trajectories, |
|
| 456 |
__global int *IdLast,int TrackPoints, |
|
| 457 |
float Mass,float InternalRadius, |
|
| 458 |
float ExternalRadius,float Angle, |
|
| 459 |
float ObserverDistance, |
|
| 460 |
int BlackBody,int AngularCamera) |
|
| 461 |
{
|
|
| 462 |
// Integer Impact Parameter ID |
|
| 463 |
int bi=get_global_id(0); |
|
| 464 |
// Integer Impact Parameter Size (half of image) |
|
| 465 |
int bmaxi=get_global_size(0); |
|
| 466 |
|
|
| 467 |
// Perform trajectory for each pixel |
|
| 468 |
|
|
| 469 |
float m,rs,ri,re,tho,ro; |
|
| 470 |
int dist,raie,q; |
|
| 471 |
|
|
| 472 |
m=Mass; |
|
| 473 |
rs=2.*m; |
|
| 474 |
ri=InternalRadius; |
|
| 475 |
re=ExternalRadius; |
|
| 476 |
ro=ObserverDistance; |
|
| 477 |
tho=Angle; |
|
| 478 |
q=-2; |
|
| 479 |
dist=AngularCamera; |
|
| 480 |
raie=BlackBody; |
|
| 481 |
|
|
| 482 |
float d,bmx,db,b,h; |
|
| 483 |
float phi,thi,phd,php,nr,r; |
|
| 484 |
int nh; |
|
| 485 |
float zp,fp; |
|
| 486 |
|
|
| 487 |
// Autosize for image |
|
| 488 |
bmx=1.25*re; |
|
| 489 |
|
|
| 490 |
// Angular step of integration |
|
| 491 |
h=4.e0f*PI/(float)TrackPoints; |
|
| 492 |
|
|
| 493 |
// impact parameter |
|
| 494 |
b=(float)bi/(float)bmaxi*bmx; |
|
| 495 |
|
|
| 496 |
float up,vp,pp,us,vs,ps; |
|
| 497 |
|
|
| 498 |
if (dist==1) |
|
| 499 |
{
|
|
| 500 |
up=0.; |
|
| 501 |
vp=1.; |
|
| 502 |
} |
|
| 503 |
else |
|
| 504 |
{
|
|
| 505 |
thi=asin(b/ro); |
|
| 506 |
up=sin(thi); |
|
| 507 |
vp=cos(thi); |
|
| 508 |
} |
|
| 509 |
|
|
| 510 |
pp=0.; |
|
| 511 |
nh=0; |
|
| 512 |
|
|
| 513 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 514 |
|
|
| 515 |
// b versus us |
|
| 516 |
float bvus=fabs(b/us); |
|
| 517 |
float bvus0=bvus; |
|
| 518 |
Trajectories[bi*TrackPoints+nh]=bvus; |
|
| 519 |
|
|
| 520 |
do |
|
| 521 |
{
|
|
| 522 |
nh++; |
|
| 523 |
pp=ps; |
|
| 524 |
up=us; |
|
| 525 |
vp=vs; |
|
| 526 |
rungekutta(&ps,&us,&vs,pp,up,vp,h,m,b); |
|
| 527 |
bvus=fabs(b/us); |
|
| 528 |
Trajectories[bi*TrackPoints+nh]=bvus; |
|
| 529 |
|
|
| 530 |
} while ((bvus>=rs)&&(bvus<=bvus0)); |
|
| 531 |
|
|
| 532 |
IdLast[bi]=nh; |
|
| 533 |
|
|
| 534 |
barrier(CLK_GLOBAL_MEM_FENCE); |
|
| 272 | 535 |
|
| 273 | 536 |
} |
| 274 | 537 |
""") |
| ... | ... | |
| 282 | 545 |
image = Image.fromarray(SigmaInt) |
| 283 | 546 |
image.save("%s.jpg" % prefix)
|
| 284 | 547 |
|
| 285 |
def BlackHole(zImage,fImage,Device):
|
|
| 548 |
def BlackHoleCL(zImage,fImage,InputCL):
|
|
| 286 | 549 |
|
| 287 | 550 |
kernel_params = {}
|
| 288 | 551 |
|
| 552 |
print(InputCL) |
|
| 553 |
|
|
| 554 |
Device=InputCL['Device'] |
|
| 555 |
GpuStyle=InputCL['GpuStyle'] |
|
| 556 |
VariableType=InputCL['VariableType'] |
|
| 557 |
Size=InputCL['Size'] |
|
| 558 |
Mass=InputCL['Mass'] |
|
| 559 |
InternalRadius=InputCL['InternalRadius'] |
|
| 560 |
ExternalRadius=InputCL['ExternalRadius'] |
|
| 561 |
Angle=InputCL['Angle'] |
|
| 562 |
ObserverDistance=InputCL['ObserverDistance'] |
|
| 563 |
Method=InputCL['Method'] |
|
| 564 |
TrackPoints=InputCL['TrackPoints'] |
|
| 565 |
|
|
| 566 |
if InputCL['BlackBody']: |
|
| 567 |
BlackBody=0 |
|
| 568 |
else: |
|
| 569 |
BlackBody=1 |
|
| 570 |
|
|
| 571 |
if InputCL['AngularCamera']: |
|
| 572 |
AngularCamera=1 |
|
| 573 |
else: |
|
| 574 |
AngularCamera=0 |
|
| 575 |
|
|
| 576 |
Trajectories=numpy.zeros((InputCL['Size']/2,InputCL['TrackPoints']), |
|
| 577 |
dtype=numpy.float32) |
|
| 578 |
IdLast=numpy.zeros(InputCL['Size']/2,dtype=numpy.int32) |
|
| 579 |
|
|
| 289 | 580 |
# Je detecte un peripherique GPU dans la liste des peripheriques |
| 290 | 581 |
Id=0 |
| 291 | 582 |
HasXPU=False |
| ... | ... | |
| 312 | 603 |
|
| 313 | 604 |
print(zImage) |
| 314 | 605 |
|
| 315 |
# Program based on Kernel2
|
|
| 606 |
TrajectoriesCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=Trajectories)
|
|
| 316 | 607 |
zImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=zImage) |
| 317 | 608 |
fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage) |
| 609 |
fImageCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=fImage) |
|
| 610 |
IdLastCL = cl.Buffer(ctx, mf.WRITE_ONLY | mf.COPY_HOST_PTR, hostbuf=IdLast) |
|
| 318 | 611 |
|
| 319 | 612 |
start_time=time.time() |
| 320 | 613 |
|
| 321 |
CLLaunch=BlackHoleCL.EachPixel(queue,(numpy.int32(zImage.shape[0]),numpy.int32(zImage.shape[1])),None,zImageCL,fImageCL) |
|
| 322 |
|
|
| 323 |
CLLaunch.wait() |
|
| 614 |
print(Trajectories.shape) |
|
| 615 |
if Method=='EachPixel': |
|
| 616 |
CLLaunch=BlackHoleCL.EachPixel(queue,(zImage.shape[0],zImage.shape[1]), |
|
| 617 |
None,zImageCL,fImageCL, |
|
| 618 |
numpy.float32(Mass), |
|
| 619 |
numpy.float32(InternalRadius), |
|
| 620 |
numpy.float32(ExternalRadius), |
|
| 621 |
numpy.float32(Angle), |
|
| 622 |
numpy.float32(ObserverDistance), |
|
| 623 |
numpy.int32(BlackBody), |
|
| 624 |
numpy.int32(AngularCamera)) |
|
| 625 |
CLLaunch.wait() |
|
| 626 |
elif Method=='TrajectoCircle': |
|
| 627 |
CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],), |
|
| 628 |
None,TrajectoriesCL,IdLastCL, |
|
| 629 |
numpy.uint32(Trajectories.shape[1]), |
|
| 630 |
numpy.float32(Mass), |
|
| 631 |
numpy.float32(InternalRadius), |
|
| 632 |
numpy.float32(ExternalRadius), |
|
| 633 |
numpy.float32(Angle), |
|
| 634 |
numpy.float32(ObserverDistance), |
|
| 635 |
numpy.int32(BlackBody), |
|
| 636 |
numpy.int32(AngularCamera)) |
|
| 637 |
|
|
| 638 |
CLLaunch=BlackHoleCL.Circle(queue,(Trajectories.shape[0], |
|
| 639 |
zImage.shape[0]*4),None, |
|
| 640 |
TrajectoriesCL,IdLastCL, |
|
| 641 |
zImageCL,fImageCL, |
|
| 642 |
numpy.uint32(Trajectories.shape[1]), |
|
| 643 |
numpy.float32(Mass), |
|
| 644 |
numpy.float32(InternalRadius), |
|
| 645 |
numpy.float32(ExternalRadius), |
|
| 646 |
numpy.float32(Angle), |
|
| 647 |
numpy.float32(ObserverDistance), |
|
| 648 |
numpy.int32(BlackBody), |
|
| 649 |
numpy.int32(AngularCamera)) |
|
| 650 |
CLLaunch.wait() |
|
| 651 |
else: |
|
| 652 |
CLLaunch=BlackHoleCL.Trajectory(queue,(Trajectories.shape[0],), |
|
| 653 |
None,TrajectoriesCL,IdLastCL, |
|
| 654 |
numpy.uint32(Trajectories.shape[1]), |
|
| 655 |
numpy.float32(Mass), |
|
| 656 |
numpy.float32(InternalRadius), |
|
| 657 |
numpy.float32(ExternalRadius), |
|
| 658 |
numpy.float32(Angle), |
|
| 659 |
numpy.float32(ObserverDistance), |
|
| 660 |
numpy.int32(BlackBody), |
|
| 661 |
numpy.int32(AngularCamera)) |
|
| 662 |
|
|
| 663 |
CLLaunch=BlackHoleCL.Pixel(queue,(zImage.shape[0], |
|
| 664 |
zImage.shape[1]),None, |
|
| 665 |
zImageCL,fImageCL,TrajectoriesCL,IdLastCL, |
|
| 666 |
numpy.uint32(Trajectories.shape[0]), |
|
| 667 |
numpy.uint32(Trajectories.shape[1]), |
|
| 668 |
numpy.float32(Mass), |
|
| 669 |
numpy.float32(InternalRadius), |
|
| 670 |
numpy.float32(ExternalRadius), |
|
| 671 |
numpy.float32(Angle), |
|
| 672 |
numpy.float32(ObserverDistance), |
|
| 673 |
numpy.int32(BlackBody), |
|
| 674 |
numpy.int32(AngularCamera)) |
|
| 675 |
CLLaunch.wait() |
|
| 676 |
|
|
| 324 | 677 |
elapsed = time.time()-start_time |
| 325 | 678 |
print("Elapsed %f: " % elapsed)
|
| 326 | 679 |
|
| 327 | 680 |
cl.enqueue_copy(queue,zImage,zImageCL).wait() |
| 328 | 681 |
cl.enqueue_copy(queue,fImage,fImageCL).wait() |
| 682 |
cl.enqueue_copy(queue,Trajectories,TrajectoriesCL).wait() |
|
| 683 |
cl.enqueue_copy(queue,IdLast,IdLastCL).wait() |
|
| 329 | 684 |
|
| 330 | 685 |
print(zImage.max()) |
| 331 | 686 |
print(fImage.max()) |
| 332 | 687 |
zImageCL.release() |
| 333 | 688 |
fImageCL.release() |
| 334 |
|
|
| 335 |
ImageOutput(zImage,"TrouNoirZ") |
|
| 336 |
ImageOutput(fImage,"TrouNoirF") |
|
| 689 |
|
|
| 690 |
TrajectoriesCL.release() |
|
| 691 |
IdLastCL.release() |
|
| 692 |
|
|
| 337 | 693 |
return(elapsed) |
| 338 | 694 |
|
| 339 | 695 |
if __name__=='__main__': |
| ... | ... | |
| 358 | 714 |
VariableType='FP32' |
| 359 | 715 |
# ? |
| 360 | 716 |
q=-2 |
| 717 |
# Method of resolution |
|
| 718 |
Method='Trajecto' |
|
| 361 | 719 |
|
| 362 |
HowToUse='%s -h [Help] -b [BlackBodyEmission] -c [AngularCamera] -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -o <ObservatorDistance> -a <AngleAboveDisc> -d <DeviceId> -g <CUDA/OpenCL> -t <FP32/FP64>' |
|
| 720 |
HowToUse='%s -h [Help] -b [BlackBodyEmission] -c [AngularCamera] -s <SizeInPixels> -m <Mass> -i <DiscInternalRadius> -x <DiscExternalRadius> -o <ObservatorDistance> -a <AngleAboveDisc> -d <DeviceId> -g <CUDA/OpenCL> -t <EachPixel/TrajecoCircle/TrajectoPixel> -v <FP32/FP64>'
|
|
| 363 | 721 |
|
| 364 | 722 |
try: |
| 365 |
opts, args = getopt.getopt(sys.argv[1:],"hbcs:m:i:x:o:a:d:g:t:",["blackbody","camera","size=","mass=","internal=","external=","observer=","angle=","device=","gpustyle=","variabletype="])
|
|
| 723 |
opts, args = getopt.getopt(sys.argv[1:],"hbcs:m:i:x:o:a:d:g:v:t:",["blackbody","camera","size=","mass=","internal=","external=","observer=","angle=","device=","gpustyle=","variabletype=","method="])
|
|
| 366 | 724 |
except getopt.GetoptError: |
| 367 | 725 |
print(HowToUse % sys.argv[0]) |
| 368 | 726 |
sys.exit(2) |
| ... | ... | |
| 409 | 767 |
Device=int(arg) |
| 410 | 768 |
elif opt in ("-g", "--gpustyle"):
|
| 411 | 769 |
GpuStyle = arg |
| 412 |
elif opt in ("-t", "--variabletype"):
|
|
| 770 |
elif opt in ("-v", "--variabletype"):
|
|
| 413 | 771 |
VariableType = arg |
| 414 | 772 |
elif opt in ("-s", "--size"):
|
| 415 | 773 |
Size = int(arg) |
| ... | ... | |
| 427 | 785 |
BlackBody = True |
| 428 | 786 |
elif opt in ("-c", "--camera"):
|
| 429 | 787 |
AngularCamera = True |
| 788 |
elif opt in ("-t", "--method"):
|
|
| 789 |
Method = arg |
|
| 430 | 790 |
|
| 431 | 791 |
print("Device Identification selected : %s" % Device)
|
| 432 | 792 |
print("GpuStyle used : %s" % GpuStyle)
|
| ... | ... | |
| 439 | 799 |
print("Angle with normal of (in radians) : %f" % Angle)
|
| 440 | 800 |
print("Black Body Disc Emission (monochromatic instead) : %s" % BlackBody)
|
| 441 | 801 |
print("Angular Camera (dimension of object instead) : %s" % AngularCamera)
|
| 802 |
print("Method of resolution : %s" % Method)
|
|
| 442 | 803 |
|
| 443 | 804 |
if GpuStyle=='CUDA': |
| 444 | 805 |
print("\nSelection of CUDA device")
|
| ... | ... | |
| 478 | 839 |
# print(Devices,Alu) |
| 479 | 840 |
|
| 480 | 841 |
# MyImage=numpy.where(numpy.random.zeros(Size,Size)>0,1,-1).astype(numpy.float32) |
| 842 |
TrackPoints=2048 |
|
| 481 | 843 |
zImage=numpy.zeros((Size,Size),dtype=numpy.float32) |
| 482 | 844 |
fImage=numpy.zeros((Size,Size),dtype=numpy.float32) |
| 483 | 845 |
|
| 484 |
print(zImage) |
|
| 485 |
|
|
| 486 |
duration=BlackHole(zImage,fImage,Device) |
|
| 846 |
InputCL={}
|
|
| 847 |
InputCL['Device']=Device |
|
| 848 |
InputCL['GpuStyle']=GpuStyle |
|
| 849 |
InputCL['VariableType']=VariableType |
|
| 850 |
InputCL['Size']=Size |
|
| 851 |
InputCL['Mass']=Mass |
|
| 852 |
InputCL['InternalRadius']=InternalRadius |
|
| 853 |
InputCL['ExternalRadius']=ExternalRadius |
|
| 854 |
InputCL['ObserverDistance']=ObserverDistance |
|
| 855 |
InputCL['Angle']=Angle |
|
| 856 |
InputCL['BlackBody']=BlackBody |
|
| 857 |
InputCL['AngularCamera']=AngularCamera |
|
| 858 |
InputCL['Method']=Method |
|
| 859 |
InputCL['TrackPoints']=TrackPoints |
|
| 860 |
|
|
| 861 |
duration=BlackHoleCL(zImage,fImage,InputCL) |
|
| 487 | 862 |
|
| 488 |
|
|
| 863 |
ImageOutput(zImage,"TrouNoirZ_%s" % Method) |
|
| 864 |
ImageOutput(fImage,"TrouNoirF_%s" % Method) |
|
Formats disponibles : Unified diff