Newer
Older
#!/usr/bin/env python3
#
# TrouNoir model using PyOpenCL or PyCUDA
#
# CC BY-NC-SA 2019 : <emmanuel.quemener@ens-lyon.fr>
#
# Part of matrix programs from: https://forge.cbp.ens-lyon.fr/svn/bench4gpu/
#
# Thanks to Andreas Klockner for PyOpenCL and PyCUDA:
# 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
# - 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 sys
import getopt
from socket import gethostname
PhysicsList = {"Einstein": 0, "Newton": 1}
return PhysicsList
#
# Blank space below to simplify debugging on OpenCL code
#
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
#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 )); // # noqa: E501
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
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,
Loading
Loading full blame...