8d4f0940
David Mayerich
ERROR in plane wa...
|
1
2
3
|
#ifndef RTS_EFIELD
#define RTS_EFIELD
|
a9275be5
David Mayerich
added vector fiel...
|
4
|
#include "../math/complex.h"
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
5
|
#include "../math/realfield.cuh"
|
a9275be5
David Mayerich
added vector fiel...
|
6
|
#include "../visualization/colormap.h"
|
a9275be5
David Mayerich
added vector fiel...
|
7
8
|
#include "../optics/planewave.h"
#include "../cuda/devices.h"
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
9
|
#include "../optics/beam.h"
|
a9275be5
David Mayerich
added vector fiel...
|
10
11
|
|
4d67ff4e
David Mayerich
started scalar fi...
|
12
13
14
|
namespace rts{
|
a9275be5
David Mayerich
added vector fiel...
|
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
|
template<typename T>
__global__ void gpu_planewave2efield(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1,
planewave<T> w, quad<T> q)
{
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;
//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;
//compute the index into the field
int i = iv*r0 + iu;
//get the current position
vec<T> p = q( (T)iu/(T)r0, (T)iv/(T)r1 );
vec<T> r(p[0], p[1], p[2]);
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
32
|
complex<T> x( 0.0f, w.k.dot(r) );
|
a9275be5
David Mayerich
added vector fiel...
|
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
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
|
if(Y == NULL) //if this is a scalar simulation
X[i] += w.E0.len() * exp(x); //use the vector magnitude as the plane wave amplitude
else
{
X[i] += w.E0[0] * exp(x);
Y[i] += w.E0[1] * exp(x);
Z[i] += w.E0[2] * exp(x);
}
}
template<typename T>
__global__ void gpu_efield_magnitude(complex<T>* X, complex<T>* Y, complex<T>* Z, unsigned int r0, unsigned int r1, T* M)
{
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;
//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;
//compute the index into the field
int i = iv*r0 + iu;
if(Y == NULL) //if this is a scalar simulation
M[i] = X[i].abs(); //compute the magnitude of the X component
else
{
M[i] = rts::vec<T>(X[i].abs(), Y[i].abs(), Z[i].abs()).len();
//M[i] = Z[i].abs();
}
}
template<typename T>
__global__ void gpu_efield_polarization(complex<T>* X, complex<T>* Y, complex<T>* Z,
unsigned int r0, unsigned int r1,
T* Px, T* Py, T* Pz)
{
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;
//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;
//compute the index into the field
int i = iv*r0 + iu;
//compute the field polarization
Px[i] = X[i].abs();
Py[i] = Y[i].abs();
Pz[i] = Z[i].abs();
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
83
84
85
86
87
88
89
90
91
|
}
/* This function computes the sum of two complex fields and stores the result in *dest
*/
template<typename T>
__global__ void gpu_efield_sum(complex<T>* dest, complex<T>* src, unsigned int r0, unsigned int r1)
{
int iu = blockIdx.x * blockDim.x + threadIdx.x;
int iv = blockIdx.y * blockDim.y + threadIdx.y;
|
a9275be5
David Mayerich
added vector fiel...
|
92
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
93
94
95
96
97
98
99
100
|
//make sure that the thread indices are in-bounds
if(iu >= r0 || iv >= r1) return;
//compute the index into the field
int i = iv*r0 + iu;
//sum the fields
dest[i] += src[i];
|
a9275be5
David Mayerich
added vector fiel...
|
101
102
|
}
|
4d67ff4e
David Mayerich
started scalar fi...
|
103
104
105
106
107
108
|
/* This class implements a discrete representation of an electromagnetic field
in 2D. The majority of this representation is done on the GPU.
*/
template<typename P>
class efield
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
109
|
protected:
|
4d67ff4e
David Mayerich
started scalar fi...
|
110
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
111
|
bool vector;
|
4d67ff4e
David Mayerich
started scalar fi...
|
112
113
114
115
116
117
118
119
120
|
//gpu pointer to the field data
rts::complex<P>* X;
rts::complex<P>* Y;
rts::complex<P>* Z;
//resolution of the discrete field
unsigned int R[2];
|
a9275be5
David Mayerich
added vector fiel...
|
121
122
123
124
125
126
127
128
129
|
//shape of the 2D field
quad<P> pos;
void from_planewave(planewave<P> p)
{
unsigned int SQRT_BLOCK = 16;
//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
130
|
|
a9275be5
David Mayerich
added vector fiel...
|
131
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
132
|
gpu_planewave2efield<P> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], p, pos);
|
a9275be5
David Mayerich
added vector fiel...
|
133
134
|
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
135
|
void init(unsigned int res0, unsigned int res1, bool _vector)
|
4d67ff4e
David Mayerich
started scalar fi...
|
136
|
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
137
|
vector = _vector; //initialize field type
|
a9275be5
David Mayerich
added vector fiel...
|
138
139
|
X = Y = Z = NULL; //initialize all pointers to NULL
|
4d67ff4e
David Mayerich
started scalar fi...
|
140
141
142
143
144
|
R[0] = res0;
R[1] = res1;
//allocate memory on the gpu
cudaMalloc(&X, sizeof(rts::complex<P>) * R[0] * R[1]);
|
a9275be5
David Mayerich
added vector fiel...
|
145
|
cudaMemset(X, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
|
4d67ff4e
David Mayerich
started scalar fi...
|
146
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
147
|
if(vector)
|
4d67ff4e
David Mayerich
started scalar fi...
|
148
|
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
149
|
std::cout<<"vector:";
|
4d67ff4e
David Mayerich
started scalar fi...
|
150
|
cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
|
a9275be5
David Mayerich
added vector fiel...
|
151
152
|
cudaMemset(Y, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
|
4d67ff4e
David Mayerich
started scalar fi...
|
153
|
cudaMalloc(&Z, sizeof(rts::complex<P>) * R[0] * R[1]);
|
a9275be5
David Mayerich
added vector fiel...
|
154
|
cudaMemset(Z, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
|
4d67ff4e
David Mayerich
started scalar fi...
|
155
156
157
|
}
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
158
159
160
161
162
163
164
165
166
167
168
169
170
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
201
|
void destroy()
{
if(X != NULL) cudaFree(X);
if(Y != NULL) cudaFree(Y);
if(Z != NULL) cudaFree(Z);
}
void shallowcpy(const rts::efield<P> & src)
{
vector = src.vector;
R[0] = src.R[0];
R[1] = src.R[1];
}
void deepcpy(const rts::efield<P> & src)
{
//perform a shallow copy
shallowcpy(src);
//allocate memory on the gpu
if(src.X != NULL)
{
cudaMalloc(&X, sizeof(rts::complex<P>) * R[0] * R[1]);
cudaMemcpy(X, src.X, sizeof(rts::complex<P>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
}
if(src.Y != NULL)
{
cudaMalloc(&Y, sizeof(rts::complex<P>) * R[0] * R[1]);
cudaMemcpy(Y, src.Y, sizeof(rts::complex<P>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
}
if(src.Z != NULL)
{
cudaMalloc(&Z, sizeof(rts::complex<P>) * R[0] * R[1]);
cudaMemcpy(Z, src.Z, sizeof(rts::complex<P>) * R[0] * R[1], cudaMemcpyDeviceToDevice);
}
}
public:
efield(unsigned int res0, unsigned int res1, bool _vector = true)
{
init(res0, res1, _vector);
pos = rts::quad<P>(rts::vec<P>(-10, 0, -10), rts::vec<P>(-10, 0, 10), rts::vec<P>(10, 0, 10));
}
|
4d67ff4e
David Mayerich
started scalar fi...
|
202
203
204
|
//destructor
~efield()
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
205
206
207
208
209
210
211
212
213
214
215
216
217
|
destroy();
}
///Clear the field - set all points to zero
void clear()
{
cudaMemset(X, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
if(vector)
{
cudaMemset(Y, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
cudaMemset(Z, 0, sizeof(rts::complex<P>) * R[0] * R[1]);
}
|
4d67ff4e
David Mayerich
started scalar fi...
|
218
219
|
}
|
a9275be5
David Mayerich
added vector fiel...
|
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
|
void position(quad<P> _p)
{
pos = _p;
}
std::string str()
{
stringstream ss;
ss<<pos<<std::endl;
ss<<"X: "<<X<<std::endl;
ss<<"Y: "<<Y<<std::endl;
ss<<"Z: "<<Z;
return ss.str();
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
236
237
238
239
240
241
242
243
|
//assignment operator: assignment from another electric field
efield<P> & operator= (const efield<P> & rhs)
{
destroy(); //destroy any previous information about this field
deepcpy(rhs); //create a deep copy
return *this; //return the current object
}
|
a9275be5
David Mayerich
added vector fiel...
|
244
245
246
247
248
249
|
//assignment operator: build an electric field from a plane wave
efield<P> & operator= (const planewave<P> & rhs)
{
clear(); //clear any previous field data
from_planewave(rhs); //create a field from the planewave
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
250
|
return *this; //return the current object
|
a9275be5
David Mayerich
added vector fiel...
|
251
252
|
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
253
254
|
//assignment operator: add an existing electric field
efield<P> & operator+= (const efield<P> & rhs)
|
a9275be5
David Mayerich
added vector fiel...
|
255
|
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
|
//if this field and the source field represent the same regions in space
if(R[0] == rhs.R[0] && R[1] == rhs.R[1] && pos == rhs.pos)
{
int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
//sum the fields
gpu_efield_sum <<<dimGrid, dimBlock>>> (X, rhs.X, R[0], R[1]);
if(Y != NULL)
{
gpu_efield_sum <<<dimGrid, dimBlock>>> (Y, rhs.Y, R[0], R[1]);
gpu_efield_sum <<<dimGrid, dimBlock>>> (Z, rhs.Z, R[0], R[1]);
}
}
else
{
std::cout<<"ERROR in efield: The two summed fields do not represent the same geometry."<<std::endl;
exit(1);
}
|
a9275be5
David Mayerich
added vector fiel...
|
279
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
280
|
return *this; //return this object
|
a9275be5
David Mayerich
added vector fiel...
|
281
282
|
}
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
283
284
285
286
287
288
289
290
291
292
293
|
efield<P> & operator= (const rts::beam<P> & rhs)
{
//get a vector of monte-carlo samples
std::vector< rts::planewave<P> > p_list = rhs.mc();
clear(); //clear any previous field data
for(unsigned int i = 0; i < p_list.size(); i++)
from_planewave(p_list[i]);
return *this;
}
|
a9275be5
David Mayerich
added vector fiel...
|
294
295
|
//return a scalar field representing field magnitude
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
296
|
realfield<P, 1, true> mag()
|
4d67ff4e
David Mayerich
started scalar fi...
|
297
|
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
298
299
|
int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
|
4d67ff4e
David Mayerich
started scalar fi...
|
300
|
|
a9275be5
David Mayerich
added vector fiel...
|
301
|
//create a scalar field to store the result
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
302
|
realfield<P, 1, true> M(R[0], R[1]);
|
a9275be5
David Mayerich
added vector fiel...
|
303
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
304
305
306
|
//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
|
a9275be5
David Mayerich
added vector fiel...
|
307
|
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
308
309
|
//compute the magnitude and store it in a scalar field
gpu_efield_magnitude<float> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], M.ptr(0));
|
a9275be5
David Mayerich
added vector fiel...
|
310
311
|
return M;
|
4d67ff4e
David Mayerich
started scalar fi...
|
312
313
|
}
|
a9275be5
David Mayerich
added vector fiel...
|
314
|
//return a vector field representing field polarization
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
315
|
realfield<P, 3, true> polarization()
|
a9275be5
David Mayerich
added vector fiel...
|
316
|
{
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
317
|
if(!vector)
|
a9275be5
David Mayerich
added vector fiel...
|
318
319
320
321
322
323
324
325
326
327
328
329
|
{
std::cout<<"ERROR: Cannot compute polarization of a scalar field."<<std::endl;
exit(1);
}
int maxThreads = rts::maxThreadsPerBlock(); //compute the optimal block size
int SQRT_BLOCK = (int)std::sqrt((float)maxThreads);
//create one thread for each detector pixel
dim3 dimBlock(SQRT_BLOCK, SQRT_BLOCK);
dim3 dimGrid((R[0] + SQRT_BLOCK -1)/SQRT_BLOCK, (R[1] + SQRT_BLOCK - 1)/SQRT_BLOCK);
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
330
|
realfield<P, 3, true> Pol(R[0], R[1]); //create a vector field to store the result
|
a9275be5
David Mayerich
added vector fiel...
|
331
332
|
//compute the polarization and store it in the vector field
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
333
|
gpu_efield_polarization<float> <<<dimGrid, dimBlock>>> (X, Y, Z, R[0], R[1], Pol.ptr(0), Pol.ptr(1), Pol.ptr(2));
|
a9275be5
David Mayerich
added vector fiel...
|
334
335
336
337
338
|
return Pol; //return the vector field
}
|
4d67ff4e
David Mayerich
started scalar fi...
|
339
340
341
342
343
|
};
|
8d4f0940
David Mayerich
ERROR in plane wa...
|
344
345
346
|
} //end namespace rts
#endif
|