forked from weft/warp
-
Notifications
You must be signed in to change notification settings - Fork 0
/
scatter_multi.cu
403 lines (338 loc) · 13 KB
/
scatter_multi.cu
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
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
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
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
#include <cuda.h>
#include <stdio.h>
#include "datadef.h"
#include "wfloat3.h"
#include "warp_device.cuh"
#include "check_cuda.h"
__global__ void scatter_multi_kernel(unsigned N, unsigned starting_index, cross_section_data* d_xsdata, particle_data* d_particles, unsigned* d_remap){
// declare shared variables
//__shared__ unsigned n_isotopes;
//__shared__ unsigned energy_grid_len;
//__shared__ unsigned total_reaction_channels;
//__shared__ unsigned* rxn_numbers;
//__shared__ unsigned* rxn_numbers_total;
//__shared__ float* energy_grid;
//__shared__ float* rxn_Q;
//__shared__ float* xs;
__shared__ float* awr;
//__shared__ float* temp;
__shared__ dist_container* dist_scatter;
__shared__ dist_container* dist_energy;
__shared__ spatial_data* space;
__shared__ unsigned* rxn;
__shared__ float* E;
__shared__ float* Q;
__shared__ unsigned* rn_bank;
//__shared__ unsigned* cellnum;
//__shared__ unsigned* matnum;
__shared__ unsigned* isonum;
//__shared__ unsigned* yield;
__shared__ float* weight;
__shared__ unsigned* index;
// have thread 0 of block copy all pointers and static info into shared memory
if (threadIdx.x == 0){
//n_isotopes = d_xsdata[0].n_isotopes;
//energy_grid_len = d_xsdata[0].energy_grid_len;
//total_reaction_channels = d_xsdata[0].total_reaction_channels;
//rxn_numbers = d_xsdata[0].rxn_numbers;
//rxn_numbers_total = d_xsdata[0].rxn_numbers_total;
//energy_grid = d_xsdata[0].energy_grid;
//rxn_Q = d_xsdata[0].Q;
//xs = d_xsdata[0].xs;
awr = d_xsdata[0].awr;
//temp = d_xsdata[0].temp;
dist_scatter = d_xsdata[0].dist_scatter;
dist_energy = d_xsdata[0].dist_energy;
space = d_particles[0].space;
rxn = d_particles[0].rxn;
E = d_particles[0].E;
Q = d_particles[0].Q;
rn_bank = d_particles[0].rn_bank;
//cellnum = d_particles[0].cellnum;
//matnum = d_particles[0].matnum;
isonum = d_particles[0].isonum;
//yield = d_particles[0].yield;
weight = d_particles[0].weight;
index = d_particles[0].index;
}
// make sure shared loads happen before anything else
__syncthreads();
// return immediately if out of bounds
int tid_in = threadIdx.x+blockIdx.x*blockDim.x;
if (tid_in >= N){return;}
//remap to active
int tid = d_remap[starting_index + tid_in];
unsigned this_rxn = rxn[ starting_index + tid_in];
// print and return if wrong
if ( this_rxn < 11 | this_rxn > 45 ){printf("multiplicity scattering kernel accessing wrong reaction @ dex %u rxn %u\n",tid, this_rxn);return;}
// check data pointers
if( dist_energy == 0x0 | dist_scatter == 0x0 ){
printf("null pointer, energy/scatter array in multiplicity scatter!,tid %u rxn %u\n",tid,this_rxn);
return;
}
//constants
//const float pi = 3.14159265359;
const float m_n = 1.00866491600; // u
//const float kb = 8.617332478e-11; // MeV/k
// load history data
wfloat3 hats_old(space[tid].xhat,space[tid].yhat,space[tid].zhat);
unsigned this_tope = isonum[ tid];
unsigned this_dex = index[ tid];
float this_E = E[ tid];
float this_weight = weight[ tid];
float this_Q = Q[ tid];
unsigned rn = rn_bank[ tid];
float this_awr = awr[ this_tope];
//float this_temp = temp[this_tope];
// pick upper or lower via stochastic mixing
dist_data this_edist, this_sdist;
dist_data sdist_lower = dist_scatter[this_dex].lower[0];
dist_data sdist_upper = dist_scatter[this_dex].upper[0];
dist_data edist_lower = dist_energy[ this_dex].lower[0];
dist_data edist_upper = dist_energy[ this_dex].upper[0];
unsigned this_law;
float f = (this_E - edist_lower.erg) / (edist_upper.erg - edist_lower.erg);
if( get_rand(&rn)>f ){
this_edist = edist_lower;
this_sdist = sdist_lower;
}
else{
this_edist = edist_upper;
this_sdist = sdist_upper;
}
this_law = this_edist.law;
// internal kernel variables
float E_target = 0;
float speed_target = sqrtf(2.0*E_target/(this_awr*m_n));
float speed_n = sqrtf(2.0*this_E/m_n);
float E_new = 0.0;
float sampled_E = 0.0;
wfloat3 v_n_cm, v_t_cm, v_n_lf, v_t_lf, v_cm, hats_new, hats_target, rotation_hat;
float mu, E0, A, R;
unsigned dist_index[1];
// ensure normalization
hats_old = hats_old / hats_old.norm2();
// make speed vectors, assume high enough energy to approximate target as stationary
v_n_lf = hats_old * speed_n;
v_t_lf = hats_target * 0.0;
// calculate v_cm
v_cm = (v_n_lf + (v_t_lf*this_awr))/(1.0+this_awr);
//transform neutron velocity into CM frame
v_n_cm = v_n_lf - v_cm;
v_t_cm = v_t_lf - v_cm;
if ( this_law == 4 ){
// sample continuous tabular
E0 = sample_continuous_tablular( this_edist.len ,
this_edist.intt ,
get_rand(&rn) ,
this_edist.var ,
this_edist.cdf,
this_edist.pdf );
//scale it to bins
sampled_E = scale_to_bins( f, E0,
this_edist.var[0], this_edist.var[ this_edist.len-1],
edist_lower.var[0], edist_lower.var[edist_lower.len-1],
edist_upper.var[0], edist_upper.var[edist_upper.len-1] );
// sample mu isotropically
mu = 2.0*get_rand(&rn)-1.0;
}
else if ( this_law == 9 ){ // evaporation spectrum
// get tabulated temperature
float t0 = edist_lower.var[0];
float t1 = edist_upper.var[0];
float U = edist_lower.cdf[0];
float e0 = edist_lower.erg;
float e1 = edist_upper.erg;
float T = 0.0;
float m = 0.0;
// interpolate T
if (e1==e0 | edist_lower.intt==1){
T = t0;
}
else if (edist_lower.intt==2){// lin-lin interpolation
T = (t1 - t0)/(e1 - e0) * (this_E - e0) + t0;
}
else{
printf("dont know what to do!\n");
}
// rejection sample
m = (this_E - U)/T;
e0 = 1.0-expf(-m);
float x = -logf(1.0-e0*get_rand(&rn)) - logf(1.0-e0*get_rand(&rn));
while ( x>m ) {
x = -logf(1.0-e0*get_rand(&rn)) - logf(1.0-e0*get_rand(&rn));
}
// mcnp5 volIII pg 2-43
sampled_E = T * x;
// isotropic mu
mu = 2.0*get_rand(&rn)-1.0;
}
else if ( this_law == 44 ){
// make sure scatter array is present
if(dist_scatter == 0x0){
printf("null pointer, scatter array in continuum !,dex %u rxn %u tope %u E %6.4E \n",this_dex,this_rxn,this_tope,this_E);
return;
}
// compute interpolation factor
if(f<0){
printf("DATA NOT WITHIN ENERGY INTERVAL tid %u rxn %u\n",tid,this_rxn);
}
// sample tabular on energy, but get index as well as value
E0 = sample_continuous_tablular( dist_index ,
this_edist.len ,
this_edist.intt ,
get_rand(&rn) ,
this_edist.var ,
this_edist.cdf,
this_edist.pdf );
//scale it to bins
sampled_E = scale_to_bins( f, E0,
this_edist.var[0], this_edist.var[ this_edist.len-1],
edist_lower.var[0], edist_lower.var[edist_lower.len-1],
edist_upper.var[0], edist_upper.var[edist_upper.len-1] );
// find correlated mu
if (this_sdist.intt==1){
A = this_sdist.var[dist_index[0]];
R = this_sdist.cdf[dist_index[0]];
}
else if (this_sdist.intt==2){
A = interpolate_linear_energy( E0,
this_edist.var[dist_index[0]],
this_edist.var[dist_index[0]+1],
this_sdist.var[dist_index[0]],
this_sdist.var[dist_index[0]+1]);
R = interpolate_linear_energy( E0,
this_edist.var[dist_index[0]],
this_edist.var[dist_index[0]+1],
this_sdist.cdf[dist_index[0]],
this_sdist.cdf[dist_index[0]+1]);
}
else{
printf("INTT=%u NOT HANDLED in law %u of multiplicity scatter!",this_sdist.law,this_sdist.intt);
}
float rn1 = get_rand(&rn);
if( get_rand(&rn)>R ){
float T = (2.0*rn1-1.0)*sinhf(A);
mu = logf(T+sqrtf(T*T+1.0))/A;
}
else{
mu = logf(rn1*expf(A)+(1.0-rn1)*expf(-A))/A;
}
}
else if ( this_law == 61 ){
// sample continuous tabular, 61 returns the the proper index depending on intt type for law 61
E0 = sample_continuous_tablular61( dist_index ,
this_edist.len ,
this_edist.intt ,
get_rand(&rn) ,
this_edist.var ,
this_edist.cdf,
this_edist.pdf );
//scale it to bins
sampled_E = scale_to_bins( f, E0,
this_edist.var[0], this_edist.var[ this_edist.len-1],
edist_lower.var[0], edist_lower.var[edist_lower.len-1],
edist_upper.var[0], edist_upper.var[edist_upper.len-1] );
// get position of data in vector and vector length
unsigned ang_position = (unsigned) this_sdist.pdf[dist_index[0]];
unsigned this_len = (unsigned) this_sdist.pdf[dist_index[0]+1] - (unsigned) this_sdist.pdf[dist_index[0]];
// sample mu from corresponding distribution
mu = sample_continuous_tablular( this_len ,
this_sdist.intt ,
get_rand(&rn) ,
&this_sdist.cdf[ ang_position ] ,
&this_sdist.cdf[ ang_position + this_sdist.len ] ,
&this_sdist.cdf[ ang_position + 2*this_sdist.len ] );
}
else if ( this_law == 66 ){
// get tabulated temperature
//unsigned nbodies = 3; // assume ONLY FOR 1002 n,2n!
unsigned this_yield1 = 2; // assume ONLY FOR 1002 n,2n!
float A = 2.9986199999999998; // assume ONLY FOR 1002 n,2n!
// if (nbodies>3) {printf("nbodies in law 66 is greater than 3!!!\n");}
// calculate things...
float Emax = (this_yield1-1.0)/this_yield1 * (A/(A+1)*this_E + this_Q);
float rn1 = get_rand(&rn);
float rn2 = get_rand(&rn);
float rn3 = get_rand(&rn);
float rn4 = get_rand(&rn);
// rejection sample random numbers
while ( (rn1*rn1 + rn2*rn2) > 1.0 ){
rn1 = get_rand(&rn);
rn2 = get_rand(&rn);
}
while ( (rn3*rn3 + rn4*rn4) > 1.0 ){
rn3 = get_rand(&rn);
rn4 = get_rand(&rn);
}
float x = -rn1 * logf(rn1*rn1 + rn2*rn2)/ (rn1*rn1 + rn2*rn2) - logf(get_rand(&rn)) ;
float y = -rn3 * logf(rn3*rn3 + rn4*rn4)/ (rn3*rn3 + rn4*rn4) - logf(get_rand(&rn)) ;
// calculate mu
sampled_E = Emax * x /( x + y );
// isotropic mu
mu = 2.0*get_rand(&rn)-1.0;
}
else{
printf("LAW %u NOT HANDLED IN MULTIPLICITY SCATTER! rxn %u\n",this_law,this_rxn);
}
// rotate direction vector
hats_old = v_n_cm / v_n_cm.norm2();
hats_old = hats_old.rotate(mu, get_rand(&rn));
// scale to sampled energy
v_n_cm = hats_old * sqrtf(2.0*sampled_E/m_n);
// transform back to L
v_n_lf = v_n_cm + v_cm;
hats_new = v_n_lf / v_n_lf.norm2();
hats_new = hats_new / hats_new.norm2(); // get higher precision, make SURE vector is length one
// calculate energy in lab frame
E_new = 0.5 * m_n * v_n_lf.dot(v_n_lf);
// check errors
if (!isfinite(sampled_E) | sampled_E < 0.0){
printf("Multiplicity scatter mis-sampled tid %i data_dex %u E %6.4E... \n",tid_in,tid,sampled_E);
}
if (!isfinite(mu) | mu < -1.0 | mu > 1.0){
printf("Multiplicity scatter mis-sampled tid %i data_dex %u mu %6.4E... \n",tid_in,tid,mu);
}
// multiply weight by multiplicity
if( this_rxn == 22 | this_rxn == 23 | this_rxn == 28 | this_rxn == 29 | this_rxn == 32 | this_rxn == 33 | this_rxn == 34 | this_rxn == 35 | this_rxn == 36 | this_rxn == 44 | this_rxn == 45){
this_weight = this_weight * 1.0;
}
else if(this_rxn == 11 | this_rxn == 16 | this_rxn == 24 | this_rxn == 30 | this_rxn == 41){
this_weight = this_weight * 2.0;
}
else if(this_rxn == 17 | this_rxn == 25 | this_rxn == 42){
this_weight = this_weight * 3.0;
}
else if(this_rxn == 37){
this_weight = this_weight * 4.0;
}
else{
printf("UNKNOWN MT NUMBER %u IN MULTIPLICITY\n",this_rxn);
}
//printf("tid %d law %u sampled_E %6.4E mu %6.4E weight %6.4E\n",tid,this_law,sampled_E,mu,this_weight);
// write results
E[ tid] = E_new;
rn_bank[tid] = rn;
weight[ tid] = this_weight;
space[tid].xhat = hats_new.x;
space[tid].yhat = hats_new.y;
space[tid].zhat = hats_new.z;
}
/**
* \brief a
* \details b
*
* @param[in] stream - CUDA stream to launch the kernel on
* @param[in] NUM_THREADS - the number of threads to run per thread block
* @param[in] N - the total number of threads to launch on the grid for multiplicity scattering
* @param[in] starting_index - starting index of the multiplicity scatter block in the remap vector
* @param[in] d_xsdata - device pointer to cross section data pointer array
* @param[in] d_particles - device pointer to particle data pointer array
* @param[in] d_remap - device pointer to data remapping vector
*/
void scatter_multi( cudaStream_t stream, unsigned NUM_THREADS, unsigned N, unsigned starting_index, cross_section_data* d_xsdata, particle_data* d_particles, unsigned* d_remap){
if(N<1){return;}
unsigned blks = ( N + NUM_THREADS - 1 ) / NUM_THREADS;
scatter_multi_kernel <<< blks, NUM_THREADS , 0 , stream >>> ( N, starting_index, d_xsdata, d_particles, d_remap );
check_cuda(cudaThreadSynchronize());
}