1 /* ----------------------------------------------------------------------
2 SPARTA - Stochastic PArallel Rarefied-gas Time-accurate Analyzer
3 http://sparta.sandia.gov
4 Steve Plimpton, sjplimp@sandia.gov, Michael Gallis, magalli@sandia.gov
5 Sandia National Laboratories
6
7 Copyright (2014) Sandia Corporation. Under the terms of Contract
8 DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
9 certain rights in this software. This software is distributed under
10 the GNU General Public License.
11
12 See the README file in the top-level SPARTA directory.
13 ------------------------------------------------------------------------- */
14
15 #ifndef SPARTA_PARTICLE_KOKKOS_H
16 #define SPARTA_PARTICLE_KOKKOS_H
17
18 #include "stdio.h"
19 #include "particle.h"
20 #include "kokkos_type.h"
21 #include "Kokkos_Random.hpp"
22 #include "rand_pool_wrap.h"
23
24 namespace SPARTA_NS {
25
26 struct struct_tdual_int_1d
27 {
28 DAT::tdual_int_1d k_view;
29 };
30
31 struct struct_tdual_float_1d
32 {
33 DAT::tdual_float_1d k_view;
34 };
35
36 struct struct_tdual_int_2d
37 {
38 DAT::tdual_int_2d k_view;
39 };
40
41 struct struct_tdual_float_2d
42 {
43 DAT::tdual_float_2d k_view;
44 };
45
46 struct TagParticleZero_cellcount{};
47 struct TagParticleCompressReactions{};
48 struct TagCopyParticleReorderDestinations{};
49 struct TagFixedMemoryReorder{};
50 struct TagFixedMemoryReorderInit{};
51 struct TagSetIcellFromPlist{};
52 struct TagParticleReorder_COPYPARTICLELIST{};
53 struct TagSetDPlistNewStyle{};
54
55 template<int NEED_ATOMICS>
56 struct TagParticleSort{};
57
58
59 class ParticleKokkos : public Particle {
60 public:
61 typedef int value_type;
62
63 // methods
64
65 ParticleKokkos(class SPARTA *);
66 ~ParticleKokkos();
67 static KOKKOS_INLINE_FUNCTION
68 int add_particle_kokkos(t_particle_1d particles, int, int, int, int,
69 double *, double *, double, double);
70 #ifndef SPARTA_KOKKOS_EXACT
71 void compress_migrate(int, int *);
72 #endif
73 void sort_kokkos();
74 void grow(int);
75 void grow_species();
76 void pre_weight() override;
77 void post_weight() override;
78 void update_class_variables();
79 int add_custom(char *, int, int);
80 void grow_custom(int, int, int);
81 void remove_custom(int);
82 void copy_custom(int, int);
83 void pack_custom(int, char *);
84 void unpack_custom(char *, int);
85
86 #ifndef SPARTA_KOKKOS_EXACT
87 typedef typename Kokkos::Random_XorShift64_Pool<DeviceType>::generator_type rand_type;
88
89 //typedef typename Kokkos::Random_XorShift1024_Pool<DeviceType>::generator_type rand_type;
90 #else
91 typedef RandWrap rand_type;
92 #endif
93
94 KOKKOS_INLINE_FUNCTION
95 double erot(int, double, rand_type &) const;
96
97 KOKKOS_INLINE_FUNCTION
98 double evib(int, double, rand_type &) const;
99
100 KOKKOS_INLINE_FUNCTION
101 void pack_custom_kokkos(int, char *) const;
102
103 KOKKOS_INLINE_FUNCTION
104 void unpack_custom_kokkos(char *, int) const;
105
106 void wrap_kokkos();
107 void sync(ExecutionSpace, unsigned int);
108 void modify(ExecutionSpace, unsigned int);
109
110 KOKKOS_INLINE_FUNCTION
111 void operator()(TagParticleZero_cellcount, const int&) const;
112
113 KOKKOS_INLINE_FUNCTION
114 void operator()(TagParticleCompressReactions, const int&) const;
115
116 template<int NEED_ATOMICS>
117 KOKKOS_INLINE_FUNCTION
118 void operator()(TagParticleSort<NEED_ATOMICS>, const int&) const;
119
120 KOKKOS_INLINE_FUNCTION
121 void operator()(TagParticleReorder_COPYPARTICLELIST, const int, int&, const bool&) const;
122
123 KOKKOS_INLINE_FUNCTION
124 void operator()(TagCopyParticleReorderDestinations, const int, int&, const bool&) const;
125
126 KOKKOS_INLINE_FUNCTION
127 void operator()(TagFixedMemoryReorder, const int&) const;
128
129 KOKKOS_INLINE_FUNCTION
130 void operator()(TagFixedMemoryReorderInit, const int&) const;
131
132 KOKKOS_INLINE_FUNCTION
133 void operator()(TagSetIcellFromPlist, const int&) const;
134
135
136 tdual_particle_1d k_particles;
137 tdual_species_1d k_species;
138 DAT::tdual_int_2d k_species2group;
139
140 typedef Kokkos::DualView<struct_tdual_int_1d*, DeviceType::array_layout, DeviceType> tdual_struct_tdual_int_1d_1d;
141 typedef Kokkos::DualView<struct_tdual_float_1d*, DeviceType::array_layout, DeviceType> tdual_struct_tdual_float_1d_1d;
142 typedef Kokkos::DualView<struct_tdual_int_2d*, DeviceType::array_layout, DeviceType> tdual_struct_tdual_int_2d_1d;
143 typedef Kokkos::DualView<struct_tdual_float_2d*, DeviceType::array_layout, DeviceType> tdual_struct_tdual_float_2d_1d;
144
145 DAT::tdual_int_1d k_ewhich,k_eicol,k_edcol;
146
147 tdual_struct_tdual_int_1d_1d k_eivec;
148 tdual_struct_tdual_float_1d_1d k_edvec;
149
150 tdual_struct_tdual_int_2d_1d k_eiarray;
151 tdual_struct_tdual_float_2d_1d k_edarray;
152
153 int sorted_kk;
154
155 inline
get_maxcellcount()156 int get_maxcellcount() {
157 return maxcellcount;
158 }
159
160 inline
set_maxcellcount(int in)161 void set_maxcellcount(int in) {
162 maxcellcount = in;
163 }
164
165 private:
166 t_particle_1d d_particles;
167 t_particle_1d d_sorted;
168 t_species_1d d_species;
169 int nParticlesWksp;
170 DAT::tdual_int_scalar k_reorder_pass;
171 DAT::t_int_scalar d_reorder_pass;
172 HAT::t_int_scalar h_reorder_pass;
173
174 int nbytes;
175 int maxcellcount,ngrid;
176 int collide_rot,vibstyle;
177 double boltz;
178
179 DAT::t_int_2d d_plist;
180 DAT::t_int_1d d_cellcount;
181
182 DAT::t_int_2d d_lists;
183 DAT::t_int_1d d_mlist;
184 DAT::t_int_1d d_slist;
185
186 HAT::t_int_2d h_lists;
187 HAT::t_int_1d h_mlist;
188 HAT::t_int_1d h_slist;
189
190 DAT::t_int_scalar d_fail_flag;
191 HAT::t_int_scalar h_fail_flag;
192
193 // work memory for reduced memory reordering
194 t_particle_1d d_pswap1;
195 t_particle_1d d_pswap2;
196 };
197
198 KOKKOS_INLINE_FUNCTION
add_particle_kokkos(t_particle_1d particles,int index,int id,int ispecies,int icell,double * x,double * v,double erot,double evib)199 int ParticleKokkos::add_particle_kokkos(t_particle_1d particles, int index, int id,
200 int ispecies, int icell, double *x, double *v, double erot, double evib)
201 {
202 OnePart tmp;
203 tmp.id = id;
204 tmp.ispecies = ispecies;
205 tmp.icell = icell;
206 tmp.x[0] = x[0];
207 tmp.x[1] = x[1];
208 tmp.x[2] = x[2];
209 tmp.v[0] = v[0];
210 tmp.v[1] = v[1];
211 tmp.v[2] = v[2];
212 tmp.erot = erot;
213 tmp.evib = evib;
214 enum{PKEEP,PINSERT,PDONE,PDISCARD,PENTRY,PEXIT,PSURF}; // same as .cpp file
215 tmp.flag = PKEEP;
216
217 int realloc = 0;
218
219 if (index < particles.extent(0)) {
220 particles[index] = tmp;
221 } else {
222 realloc = 1;
223 }
224
225 return realloc;
226 }
227
228
229 /* ----------------------------------------------------------------------
230 generate random rotational energy for a particle
231 only a function of species index and species properties
232 ------------------------------------------------------------------------- */
233
234 KOKKOS_INLINE_FUNCTION
erot(int isp,double temp_thermal,rand_type & erandom)235 double ParticleKokkos::erot(int isp, double temp_thermal, rand_type &erandom) const
236 {
237 double eng,a,erm,b;
238
239 if (!collide_rot) return 0.0;
240 if (d_species[isp].rotdof < 2) return 0.0;
241
242 if (d_species[isp].rotdof == 2)
243 eng = -log(erandom.drand()) * boltz * temp_thermal;
244 else {
245 a = 0.5*d_species[isp].rotdof-1.0;
246 while (1) {
247 // energy cut-off at 10 kT
248 erm = 10.0*erandom.drand();
249 b = pow(erm/a,a) * exp(a-erm);
250 if (b > erandom.drand()) break;
251 }
252 eng = erm * boltz * temp_thermal;
253 }
254
255 return eng;
256 }
257
258 /* ----------------------------------------------------------------------
259 generate random vibrational energy for a particle
260 only a function of species index and species properties
261 ------------------------------------------------------------------------- */
262
263 KOKKOS_INLINE_FUNCTION
evib(int isp,double temp_thermal,rand_type & erandom)264 double ParticleKokkos::evib(int isp, double temp_thermal, rand_type &erandom) const
265 {
266 double eng,a,erm,b;
267
268 enum{NONE,DISCRETE,SMOOTH}; // several files
269 if (vibstyle == NONE || d_species[isp].vibdof < 2) return 0.0;
270
271 eng = 0.0;
272 if (vibstyle == DISCRETE && d_species[isp].vibdof == 2) {
273 int ivib = static_cast<int> (-log(erandom.drand()) * temp_thermal /
274 d_species[isp].vibtemp[0]);
275 eng = ivib * boltz * d_species[isp].vibtemp[0];
276 } else if (vibstyle == SMOOTH || d_species[isp].vibdof >= 2) {
277 if (d_species[isp].vibdof == 2)
278 eng = -log(erandom.drand()) * boltz * temp_thermal;
279 else if (d_species[isp].vibdof > 2) {
280 a = 0.5*d_species[isp].vibdof-1.;
281 while (1) {
282 // energy cut-off at 10 kT
283 erm = 10.0*erandom.drand();
284 b = pow(erm/a,a) * exp(a-erm);
285 if (b > erandom.drand()) break;
286 }
287 eng = erm * boltz * temp_thermal;
288 }
289 }
290
291 return eng;
292 }
293
294 KOKKOS_INLINE_FUNCTION
pack_custom_kokkos(int n,char * buf)295 void ParticleKokkos::pack_custom_kokkos(int n, char *buf) const
296 {
297 int i,j;
298 char *ptr = buf;
299
300 if (ncustom_ivec) {
301 for (i = 0; i < ncustom_ivec; i++) {
302 memcpy(ptr,&(k_eivec.d_view(i).k_view.d_view(n)),sizeof(int));
303 ptr += sizeof(int);
304 }
305 }
306 if (ncustom_iarray) {
307 for (i = 0; i < ncustom_iarray; i++) {
308 const int ncols = k_eicol.d_view[i];
309 for (j = 0; j < ncols; j++) {
310 memcpy(ptr,&(k_eiarray.d_view(i).k_view.d_view(n,j)),sizeof(int));
311 ptr += sizeof(int);
312 }
313 }
314 }
315
316 ptr = ROUNDUP(ptr);
317
318 if (ncustom_dvec) {
319 for (i = 0; i < ncustom_dvec; i++) {
320 memcpy(ptr,&(k_edvec.d_view(i).k_view.d_view(n)),sizeof(double));
321 ptr += sizeof(double);
322 }
323 }
324 if (ncustom_darray) {
325 for (i = 0; i < ncustom_darray; i++) {
326 const int ncols = k_edcol.d_view[i];
327 for (j = 0; j < ncols; j++) {
328 memcpy(ptr,&(k_edarray.d_view(i).k_view.d_view(n,j)),sizeof(double));
329 ptr += sizeof(double);
330 }
331 }
332 }
333 }
334
335 KOKKOS_INLINE_FUNCTION
unpack_custom_kokkos(char * buf,int n)336 void ParticleKokkos::unpack_custom_kokkos(char *buf, int n) const
337 {
338 int i,j;
339 char *ptr = buf;
340
341 if (ncustom_ivec) {
342 for (i = 0; i < ncustom_ivec; i++) {
343 memcpy(&(k_eivec.d_view(i).k_view.d_view(n)),ptr,sizeof(int));
344 ptr += sizeof(int);
345 }
346 }
347 if (ncustom_iarray) {
348 for (i = 0; i < ncustom_iarray; i++) {
349 const int ncols = k_eicol.d_view[i];
350 for (j = 0; j < ncols; j++) {
351 memcpy(&(k_eiarray.d_view(i).k_view.d_view(n,j)),ptr,sizeof(int));
352 ptr += sizeof(int);
353 }
354 }
355 }
356
357 ptr = ROUNDUP(ptr);
358
359 if (ncustom_dvec) {
360 for (i = 0; i < ncustom_dvec; i++) {
361 memcpy(&(k_edvec.d_view(i).k_view.d_view(n)),ptr,sizeof(double));
362 ptr += sizeof(double);
363 }
364 }
365 if (ncustom_darray) {
366 for (i = 0; i < ncustom_darray; i++) {
367 const int ncols = k_edcol.d_view[i];
368 for (j = 0; j < ncols; j++) {
369 memcpy(&(k_edarray.d_view(i).k_view.d_view(n,j)),ptr,sizeof(double));
370 ptr += sizeof(double);
371 }
372 }
373 }
374 }
375
376
377 }
378
379 #endif
380
381 /* ERROR/WARNING messages:
382
383 */
384