0c9bf8ae
dmayerich
Case-sensitive er...
|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
|
__device__ double g(double v0, double v1)
{
return (v0 + v1)*log(abs(v0+v1)) + (v0-v1)*log(abs(v0-v1));
}
__device__ double hfin(double v0, double v1, double dv)
{
double e = 0.001;
double t0 = g(v0+e, v1-dv)/dv;
double t1 = 2*g(v0+e, v1)/dv;
double t2 = g(v0+e, v1+dv)/dv;
return -1.0/PI * (t0 - t1 + t2);
}
__global__ void devKramersKronig(double* gpuN, double* gpuK, int numVals, double nuStart, double nuEnd, double nOffset)
{
int i = blockIdx.x * blockDim.x + threadIdx.x;
if(i >= numVals) return;
double nuDelta = (nuEnd - nuStart)/(numVals - 1);
double nu = nuStart + i*nuDelta;
double n = 0.0;
double jNu;
double jK;
for(int j=1; j<numVals-1; j++)
{
jNu = nuStart + j*nuDelta;
jK = gpuK[j];
n += hfin(nu, jNu, nuDelta) * jK;
}
gpuN[i] = n + nOffset;
}
void cudaKramersKronig(double* cpuN, double* cpuK, int nVals, double nuStart, double nuEnd, double nOffset)
{
double* gpuK;
HANDLE_ERROR(cudaMalloc(&gpuK, sizeof(double)*nVals));
HANDLE_ERROR(cudaMemcpy(gpuK, cpuK, sizeof(double)*nVals, cudaMemcpyHostToDevice));
double* gpuN;
HANDLE_ERROR(cudaMalloc(&gpuN, sizeof(double)*nVals));
dim3 block(BLOCK_SIZE*BLOCK_SIZE);
dim3 grid(nVals/block.x + 1);
devKramersKronig<<<grid, block>>>(gpuN, gpuK, nVals, nuStart, nuEnd, nOffset);
HANDLE_ERROR(cudaMemcpy(cpuN, gpuN, sizeof(double)*nVals, cudaMemcpyDeviceToHost));
//free resources
HANDLE_ERROR(cudaFree(gpuK));
HANDLE_ERROR(cudaFree(gpuN));
}
|
39a7d6e9
dmayerich
Added dialog supp...
|
57
58
|
__global__ void devComputeSpectrum(double* I, double2* B, double* alpha, int Nl,
int nSamples, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO)
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
59
|
{
|
39a7d6e9
dmayerich
Added dialog supp...
|
60
|
int i = blockIdx.x * blockDim.x + threadIdx.x;
|
4198e3be
dmayerich
Added polystyrene.
|
61
|
if(i >= nLambda)
|
39a7d6e9
dmayerich
Added dialog supp...
|
62
|
return;
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
63
64
65
66
67
|
//compute the delta-theta value
double dTheta = (oThetaO - oThetaI)/nSamples;
//allocate space for the Legendre polynomials
|
39a7d6e9
dmayerich
Added dialog supp...
|
68
|
double Ptheta[2];
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
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
|
double cosTheta, theta;
cuDoubleComplex Us;
cuDoubleComplex UsSample;
cuDoubleComplex U;
//cuComplex Ui;
//Ui.x = 2*PI;
//Ui.y = 0.0;
cuDoubleComplex numer;
numer.x = 0.0;
cuDoubleComplex exp_numer;
cuDoubleComplex iL;
cuDoubleComplex imag;
imag.x = 0.0; imag.y = 1.0;
double realFac;
cuDoubleComplex complexFac;
double PlTheta;
double Isum = 0.0;
//float maxVal = 0;
//float val;
for(int iTheta = 0; iTheta < nSamples; iTheta++)
{
//calculate theta
theta = iTheta * dTheta + oThetaI;
cosTheta = cos(theta);
//initialize the theta Legendre polynomial
Ptheta[0] = 1.0;
Ptheta[1] = cosTheta;
//initialize the scattered field
Us.x = Us.y = 0.0;
iL.x = 1.0;
iL.y = 0.0;
for(int l = 0; l<Nl; l++)
{
//compute the theta legendre polynomial
if(l == 0)
PlTheta = Ptheta[0];
else if(l == 1)
PlTheta = Ptheta[1];
else
{
PlTheta = ((2*l - 1)*cosTheta*Ptheta[1] - (l - 1)*Ptheta[0])/l;
Ptheta[0] = Ptheta[1];
Ptheta[1] = PlTheta;
}
//compute the real components of the scattered field
realFac = alpha[l] * PlTheta;
//compute the complex components of the scattered field
numer.x = 0.0;
numer.y = -(l*PI)/2.0;
exp_numer = cExp(numer);
complexFac = cMult(B[Nl * i + l], exp_numer);
complexFac = cMult(complexFac, iL);
|
39a7d6e9
dmayerich
Added dialog supp...
|
127
|
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
128
129
130
131
132
133
134
135
|
//combine the real and complex components
UsSample = cMult(complexFac, realFac);
Us = cAdd(Us, UsSample);
//increment the imaginary exponent i^l
iL = cMult(iL, imag);
|
39a7d6e9
dmayerich
Added dialog supp...
|
136
|
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
|
}
//sum the scattered and incident fields
if(theta >= cThetaI && theta <= cThetaO)
U = cAdd(Us, 2*PI);
else
U = Us;
Isum += (U.x*U.x + U.y*U.y) * sin(theta) * 2 * PI * dTheta;
}
I[i] = Isum;
}
void cudaComputeSpectrum(double* cpuI, double* cpuB, double* cpuAlpha,
int Nl, int nLambda, double oThetaI, double oThetaO, double cThetaI, double cThetaO, int nSamples)
{
//copy everything to the GPU
double2* gpuB;
HANDLE_ERROR(cudaMalloc(&gpuB, sizeof(double2) * nLambda * Nl));
HANDLE_ERROR(cudaMemcpy(gpuB, cpuB, sizeof(double2) * nLambda * Nl, cudaMemcpyHostToDevice));
|
4198e3be
dmayerich
Added polystyrene.
|
158
|
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
159
160
161
162
163
|
double* gpuAlpha;
HANDLE_ERROR(cudaMalloc(&gpuAlpha, sizeof(double) * Nl));
HANDLE_ERROR(cudaMemcpy(gpuAlpha, cpuAlpha, sizeof(double) * Nl, cudaMemcpyHostToDevice));
double* gpuI;
|
39a7d6e9
dmayerich
Added dialog supp...
|
164
165
166
|
HANDLE_ERROR(cudaMalloc(&gpuI, sizeof(double) * nLambda));
HANDLE_ERROR(cudaMemset(gpuI, 0, sizeof(double) * nLambda));
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
167
168
169
170
171
172
173
|
//call the kernel to compute the spectrum
dim3 block(BLOCK_SIZE*BLOCK_SIZE);
dim3 grid(nLambda/block.x + 1);
//devComputeSpectrum
devComputeSpectrum<<<grid, block>>>(gpuI, (double2*)gpuB, gpuAlpha, Nl,
|
39a7d6e9
dmayerich
Added dialog supp...
|
174
|
nSamples, nLambda, oThetaI, oThetaO, cThetaI, cThetaO);
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
175
176
177
|
HANDLE_ERROR(cudaMemcpy(cpuI, gpuI, sizeof(double) * nLambda, cudaMemcpyDeviceToHost));
|
39a7d6e9
dmayerich
Added dialog supp...
|
178
179
|
//printf("Final array value: %f\n", cpuI[nLambda-1]);
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
180
181
182
183
|
HANDLE_ERROR(cudaFree(gpuB));
HANDLE_ERROR(cudaFree(gpuAlpha));
HANDLE_ERROR(cudaFree(gpuI));
|
0c9bf8ae
dmayerich
Case-sensitive er...
|
184
185
|
|
39a7d6e9
dmayerich
Added dialog supp...
|
186
187
|
}
|