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 #include "spatype.h"
16 #include "mpi.h"
17 #include "math.h"
18 #include "stdlib.h"
19 #include "string.h"
20 #include "update_kokkos.h"
21 #include "math_const.h"
22 #include "particle_kokkos.h"
23 #include "modify.h"
24 #include "fix.h"
25 #include "compute.h"
26 #include "domain.h"
27 #include "comm_kokkos.h"
28 #include "collide.h"
29 #include "collide_vss_kokkos.h"
30 #include "grid_kokkos.h"
31 #include "surf_kokkos.h"
32 #include "surf_collide.h"
33 #include "surf_react.h"
34 #include "output.h"
35 #include "geometry_kokkos.h"
36 #include "random_mars.h"
37 #include "timer.h"
38 #include "math_extra.h"
39 #include "memory_kokkos.h"
40 #include "error.h"
41 #include <unistd.h>
42 #include "kokkos.h"
43 #include "sparta_masks.h"
44 #include "surf_collide_specular_kokkos.h"
45
46 using namespace SPARTA_NS;
47
48 enum{XLO,XHI,YLO,YHI,ZLO,ZHI,INTERIOR}; // same as Domain
49 enum{PERIODIC,OUTFLOW,REFLECT,SURFACE,AXISYM}; // same as Domain
50 //enum{OUTSIDE,INSIDE,ONSURF2OUT,ONSURF2IN}; // several files
51 enum{PKEEP,PINSERT,PDONE,PDISCARD,PENTRY,PEXIT,PSURF}; // several files
52 enum{NCHILD,NPARENT,NUNKNOWN,NPBCHILD,NPBPARENT,NPBUNKNOWN,NBOUND}; // Grid
53 enum{TALLYAUTO,TALLYREDUCE,TALLYLOCAL}; // same as Surf
54
55 #define MAXSTUCK 20
56 #define EPSPARAM 1.0e-7
57
58 // either set ID or PROC/INDEX, set other to -1
59
60 //#define MOVE_DEBUG 1 // un-comment to debug one particle
61 #define MOVE_DEBUG_ID 308143534 // particle ID
62 #define MOVE_DEBUG_PROC -1 // owning proc
63 #define MOVE_DEBUG_INDEX -1 // particle index on owning proc
64 #define MOVE_DEBUG_STEP 4107 // timestep
65
66 #define VAL_1(X) X
67 #define VAL_2(X) VAL_1(X), VAL_1(X)
68
69 /* ---------------------------------------------------------------------- */
70
UpdateKokkos(SPARTA * sparta)71 UpdateKokkos::UpdateKokkos(SPARTA *sparta) : Update(sparta),
72 grid_kk_copy(sparta),
73 domain_kk_copy(sparta),
74 //// Virtual functions are not yet supported on the GPU, which leads to pain:
75 sc_kk_specular_copy{VAL_2(KKCopy<SurfCollideSpecularKokkos>(sparta))},
76 sc_kk_diffuse_copy{VAL_2(KKCopy<SurfCollideDiffuseKokkos>(sparta))},
77 sc_kk_vanish_copy{VAL_2(KKCopy<SurfCollideVanishKokkos>(sparta))},
78 sc_kk_piston_copy{VAL_2(KKCopy<SurfCollidePistonKokkos>(sparta))},
79 sc_kk_transparent_copy{VAL_2(KKCopy<SurfCollideTransparentKokkos>(sparta))},
80 blist_active_copy{VAL_2(KKCopy<ComputeBoundaryKokkos>(sparta))},
81 slist_active_copy{VAL_2(KKCopy<ComputeSurfKokkos>(sparta))}
82 {
83
84 // use 1D view for scalars to reduce GPU memory operations
85
86 d_scalars = t_int_12("collide:scalars");
87 h_scalars = t_host_int_12("collide:scalars_mirror");
88
89 d_ncomm_one = Kokkos::subview(d_scalars,0);
90 d_nexit_one = Kokkos::subview(d_scalars,1);
91 d_nboundary_one = Kokkos::subview(d_scalars,2);
92 d_nmigrate = Kokkos::subview(d_scalars,3);
93 d_entryexit = Kokkos::subview(d_scalars,4);
94 d_ntouch_one = Kokkos::subview(d_scalars,5);
95 d_nscheck_one = Kokkos::subview(d_scalars,6);
96 d_nscollide_one = Kokkos::subview(d_scalars,7);
97 d_nreact_one = Kokkos::subview(d_scalars,8);
98 d_nstuck = Kokkos::subview(d_scalars,9);
99 d_naxibad = Kokkos::subview(d_scalars,10);
100 d_error_flag = Kokkos::subview(d_scalars,11);
101
102 h_ncomm_one = Kokkos::subview(h_scalars,0);
103 h_nexit_one = Kokkos::subview(h_scalars,1);
104 h_nboundary_one = Kokkos::subview(h_scalars,2);
105 h_nmigrate = Kokkos::subview(h_scalars,3);
106 h_entryexit = Kokkos::subview(h_scalars,4);
107 h_ntouch_one = Kokkos::subview(h_scalars,5);
108 h_nscheck_one = Kokkos::subview(h_scalars,6);
109 h_nscollide_one = Kokkos::subview(h_scalars,7);
110 h_nreact_one = Kokkos::subview(h_scalars,8);
111 h_nstuck = Kokkos::subview(h_scalars,9);
112 h_naxibad = Kokkos::subview(h_scalars,10);
113 h_error_flag = Kokkos::subview(h_scalars,11);
114
115 nboundary_tally = 0;
116 }
117
118 /* ---------------------------------------------------------------------- */
119
~UpdateKokkos()120 UpdateKokkos::~UpdateKokkos()
121 {
122 if (copymode) return;
123
124 memoryKK->destroy_kokkos(k_mlist,mlist);
125 mlist = NULL;
126
127 grid_kk_copy.uncopy();
128 domain_kk_copy.uncopy();
129
130 for (int i=0; i<KOKKOS_MAX_SURF_COLL_PER_TYPE; i++) {
131 sc_kk_specular_copy[i].uncopy();
132 sc_kk_diffuse_copy[i].uncopy();
133 sc_kk_vanish_copy[i].uncopy();
134 sc_kk_piston_copy[i].uncopy();
135 sc_kk_transparent_copy[i].uncopy();
136 }
137
138 for (int i=0; i<KOKKOS_MAX_BLIST; i++) {
139 blist_active_copy[i].uncopy();
140 }
141
142 for (int i=0; i<KOKKOS_MAX_SLIST; i++) {
143 slist_active_copy[i].uncopy();
144 }
145 }
146
147 /* ---------------------------------------------------------------------- */
148
init()149 void UpdateKokkos::init()
150 {
151 // init the UpdateKokkos class if performing a run, else just return
152 // only set first_update if a run is being performed
153
154 if (runflag == 0) return;
155 first_update = 1;
156
157 // choose the appropriate move method
158
159 moveptr = NULL;
160 if (domain->dimension == 3) {
161 if (surf->exist) moveptr = &UpdateKokkos::move<3,1>;
162 else moveptr = &UpdateKokkos::move<3,0>;
163 } else if (domain->axisymmetric) {
164 if (surf->exist) moveptr = &UpdateKokkos::move<1,1>;
165 else moveptr = &UpdateKokkos::move<1,0>;
166 } else if (domain->dimension == 2) {
167 if (surf->exist) moveptr = &UpdateKokkos::move<2,1>;
168 else moveptr = &UpdateKokkos::move<2,0>;
169 }
170
171 // check gravity vector
172
173 if (domain->dimension == 2 && gravity[2] != 0.0)
174 error->all(FLERR,"Gravity in z not allowed for 2d");
175 if (domain->axisymmetric && gravity[1] != 0.0)
176 error->all(FLERR,"Gravity in y not allowed for axi-symmetric model");
177
178 // moveperturb method is set if particle motion is perturbed
179
180 gravity_3d_flag = gravity_2d_flag = 0;
181 moveperturb = NULL;
182 if (gravity[0] != 0.0 || gravity[1] != 0.0 || gravity[2] != 0.0) {
183 if (domain->dimension == 3) {
184 moveperturb = &UpdateKokkos::gravity3d;
185 gravity_3d_flag = 1;
186 }
187 if (domain->dimension == 2){
188 moveperturb = &UpdateKokkos::gravity2d;
189 gravity_2d_flag = 1;
190 }
191 }
192
193 if (moveperturb) perturbflag = 1;
194 else perturbflag = 0;
195
196 }
197
198 /* ---------------------------------------------------------------------- */
199
setup()200 void UpdateKokkos::setup()
201 {
202 ParticleKokkos* particle_kk = (ParticleKokkos*) particle;
203 GridKokkos* grid_kk = (GridKokkos*) grid;
204 SurfKokkos* surf_kk = (SurfKokkos*) surf;
205
206 particle_kk->sync(Device,ALL_MASK);
207
208 if (sparta->kokkos->prewrap) {
209
210 // particle
211
212 particle_kk->wrap_kokkos();
213
214 // grid
215
216 grid_kk->wrap_kokkos();
217 grid_kk->update_hash();
218
219 // surf
220
221 if (surf->exist) {
222 surf_kk->wrap_kokkos();
223 }
224
225 sparta->kokkos->prewrap = 0;
226 } else {
227 grid_kk->modify(Host,ALL_MASK);
228 grid_kk->update_hash();
229 if (surf->exist) {
230 surf_kk->modify(Host,ALL_MASK);
231 grid_kk->wrap_kokkos_graphs();
232 }
233 }
234
235 Update::setup(); // must come after prewrap since computes are called by setup()
236
237 //// For MPI debugging
238 //
239 // volatile int i = 0;
240 // char hostname[256];
241 // gethostname(hostname, sizeof(hostname));
242 // printf("PID %d on %s ready for attach, i = %i\n", getpid(), hostname, i);
243 // fflush(stdout);
244 // sleep(30);
245 // printf("Continuing...\n");
246 }
247
248 /* ---------------------------------------------------------------------- */
249
run(int nsteps)250 void UpdateKokkos::run(int nsteps)
251 {
252 sparta->kokkos->auto_sync = 0;
253
254 ParticleKokkos* particle_kk = (ParticleKokkos*) particle;
255
256 int n_start_of_step = modify->n_start_of_step;
257 int n_end_of_step = modify->n_end_of_step;
258
259 // cellweightflag = 1 if grid-based particle weighting is ON
260
261 int cellweightflag = 0;
262 if (grid->cellweightflag) cellweightflag = 1;
263
264 // loop over timesteps
265
266 for (int i = 0; i < nsteps; i++) {
267
268 ntimestep++;
269 if (collide_react) collide_react_update();
270 if (bounce_tally) bounce_set(ntimestep);
271
272 timer->stamp();
273
274 // dynamic parameter updates
275
276 if (dynamic) dynamic_update();
277
278 // start of step fixes
279
280 if (n_start_of_step) {
281 modify->start_of_step();
282 timer->stamp(TIME_MODIFY);
283 }
284
285 // move particles
286
287 if (cellweightflag) particle->pre_weight();
288 (this->*moveptr)();
289 timer->stamp(TIME_MOVE);
290
291 // communicate particles
292
293 if (nmigrate) {
294 k_mlist_small = Kokkos::subview(k_mlist,std::make_pair(0,nmigrate));
295 k_mlist_small.sync_host();
296 }
297 auto mlist_small = k_mlist_small.h_view.data();
298
299 ((CommKokkos*)comm)->migrate_particles(nmigrate,mlist_small,k_mlist_small.d_view);
300 if (cellweightflag) particle->post_weight();
301 timer->stamp(TIME_COMM);
302
303 if (collide) {
304 particle_kk->sort_kokkos();
305 timer->stamp(TIME_SORT);
306
307 //CollideVSSKokkos* collide_kk = (CollideVSSKokkos*) collide;
308 //collide_kk->wrap_kokkos_sort();
309
310 collide->collisions();
311 timer->stamp(TIME_COLLIDE);
312 }
313
314 // diagnostic fixes
315
316 if (n_end_of_step) {
317 modify->end_of_step();
318 timer->stamp(TIME_MODIFY);
319 }
320
321 // all output
322
323 if (ntimestep == output->next) {
324 particle_kk->sync(Host,ALL_MASK);
325 output->write(ntimestep);
326 timer->stamp(TIME_OUTPUT);
327 }
328 }
329 sparta->kokkos->auto_sync = 1;
330
331 particle_kk->sync(Host,ALL_MASK);
332
333 }
334
335 //make randomread versions of d_particles?
336
337 /* ----------------------------------------------------------------------
338 advect particles thru grid
339 DIM = 2/3 for 2d/3d, 1 for 2d axisymmetric
340 SURF = 0/1 for no surfs or surfs
341 use multiple iterations of move/comm if necessary
342 ------------------------------------------------------------------------- */
343
move()344 template < int DIM, int SURF > void UpdateKokkos::move()
345 {
346 //bool hitflag;
347 //int m,icell,icell_original,nmask,outface,bflag,nflag,pflag,itmp;
348 //int side,minside,minsurf,nsurf,cflag,isurf,exclude,stuck_iterate;
349 //int pstart,pstop,entryexit,any_entryexit;
350 //cellint *neigh;
351 //double dtremain,frac,newfrac,param,minparam,rnew,dtsurf,tc,tmp;
352 //double xnew[3],xhold[3],xc[3],vc[3],minxc[3],minvc[3];
353 //double *x,*v,*lo,*hi;
354 //Particle::OnePart iorig;
355 //Particle::OnePart *particles;
356 //Particle::OnePart *ipart,*jpart;
357
358 int pstart,pstop,entryexit,any_entryexit;
359
360 // extend migration list if necessary
361
362 int nlocal = particle->nlocal;
363 int maxlocal = particle->maxlocal;
364
365 if (nlocal > maxmigrate) {
366 maxmigrate = maxlocal;
367 memoryKK->destroy_kokkos(k_mlist,mlist);
368 memoryKK->create_kokkos(k_mlist,mlist,maxmigrate,"particle:mlist");
369 }
370
371 // counters
372
373 niterate = 0;
374 ntouch_one = ncomm_one = 0;
375 nboundary_one = nexit_one = 0;
376 nscheck_one = nscollide_one = 0;
377 surf->nreact_one = 0;
378
379 if (!sparta->kokkos->need_atomics || sparta->kokkos->atomic_reduction) {
380 h_ntouch_one() = 0;
381 h_nexit_one() = 0;
382 h_nboundary_one() = 0;
383 h_ncomm_one() = 0;
384 h_nscheck_one() = 0;
385 h_nscollide_one() = 0;
386 h_nreact_one() = 0;
387 }
388
389 h_error_flag() = 0;
390
391 // move/migrate iterations
392
393 dt = update->dt;
394
395 ParticleKokkos* particle_kk = ((ParticleKokkos*)particle);
396
397 while (1) {
398
399 // loop over particles
400 // first iteration = all my particles
401 // subsequent iterations = received particles
402
403 niterate++;
404
405 d_particles = particle_kk->k_particles.d_view;
406
407 GridKokkos* grid_kk = ((GridKokkos*)grid);
408 d_cells = grid_kk->k_cells.d_view;
409 d_sinfo = grid_kk->k_sinfo.d_view;
410 d_pcells = grid_kk->k_pcells.d_view;
411
412 d_csurfs = grid_kk->d_csurfs;
413 d_csplits = grid_kk->d_csplits;
414 d_csubs = grid_kk->d_csubs;
415
416 if (surf->exist) {
417 SurfKokkos* surf_kk = ((SurfKokkos*)surf);
418 surf_kk->sync(Device,ALL_MASK);
419 d_lines = surf_kk->k_lines.d_view;
420 d_tris = surf_kk->k_tris.d_view;
421 }
422
423 particle_kk->sync(Device,PARTICLE_MASK);
424 grid_kk->sync(Device,CELL_MASK|PCELL_MASK|SINFO_MASK|PLEVEL_MASK);
425
426 // may be able to move this outside of the while loop
427 grid_kk_copy.copy(grid_kk);
428 domain_kk_copy.copy((DomainKokkos*)domain);
429
430 if (surf->nsc > KOKKOS_TOT_SURF_COLL)
431 error->all(FLERR,"Kokkos currently supports two instances of each surface collide method");
432
433 if (surf->nsc > 0) {
434 int nspec,ndiff,nvan,npist,ntrans;
435 nspec = ndiff = nvan = npist = ntrans = 0;
436 for (int n = 0; n < surf->nsc; n++) {
437 if (strcmp(surf->sc[n]->style,"specular") == 0) {
438 sc_kk_specular_copy[nspec].copy((SurfCollideSpecularKokkos*)(surf->sc[n]));
439 sc_type_list[n] = 0;
440 sc_map[n] = nspec;
441 nspec++;
442 } else if (strcmp(surf->sc[n]->style,"diffuse") == 0) {
443 sc_kk_diffuse_copy[ndiff].copy((SurfCollideDiffuseKokkos*)(surf->sc[n]));
444 sc_kk_diffuse_copy[ndiff].obj.pre_collide();
445 sc_type_list[n] = 1;
446 sc_map[n] = ndiff;
447 ndiff++;
448 } else if (strcmp(surf->sc[n]->style,"vanish") == 0) {
449 sc_kk_vanish_copy[nvan].copy((SurfCollideVanishKokkos*)(surf->sc[n]));
450 sc_type_list[n] = 2;
451 sc_map[n] = nvan;
452 nvan++;
453 } else if (strcmp(surf->sc[n]->style,"piston") == 0) {
454 sc_kk_piston_copy[npist].copy((SurfCollidePistonKokkos*)(surf->sc[n]));
455 sc_type_list[n] = 3;
456 sc_map[n] = npist;
457 npist++;
458 } else if (strcmp(surf->sc[n]->style,"transparent") == 0) {
459 sc_kk_transparent_copy[ntrans].copy((SurfCollideTransparentKokkos*)(surf->sc[n]));
460 sc_type_list[n] = 4;
461 sc_map[n] = ntrans;
462 ntrans++;
463 } else {
464 error->all(FLERR,"Unknown Kokkos surface collide method");
465 }
466 }
467 if (nspec > KOKKOS_MAX_SURF_COLL_PER_TYPE || ndiff > KOKKOS_MAX_SURF_COLL_PER_TYPE ||
468 nvan > KOKKOS_MAX_SURF_COLL_PER_TYPE || npist > KOKKOS_MAX_SURF_COLL_PER_TYPE ||
469 ntrans > KOKKOS_MAX_SURF_COLL_PER_TYPE)
470 error->all(FLERR,"Kokkos currently supports two instances of each surface collide method");
471 }
472
473 if (surf->nsr)
474 error->all(FLERR,"Cannot (yet) use surface reactions with Kokkos");
475
476 h_nmigrate() = 0;
477 h_entryexit() = 0;
478
479 nmigrate = 0;
480 entryexit = 0;
481
482 if (niterate == 1) {
483 pstart = 0;
484 pstop = nlocal;
485 }
486
487 UPDATE_REDUCE reduce;
488
489 /* ATOMIC_REDUCTION: 1 = use atomics
490 0 = don't need atomics
491 -1 = use parallel_reduce
492 */
493
494 Kokkos::deep_copy(d_scalars,h_scalars);
495
496 //k_mlist.sync_device();
497 copymode = 1;
498 if (!sparta->kokkos->need_atomics)
499 Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagUpdateMove<DIM,SURF,0> >(pstart,pstop),*this);
500 else if (sparta->kokkos->atomic_reduction)
501 Kokkos::parallel_for(Kokkos::RangePolicy<DeviceType, TagUpdateMove<DIM,SURF,1> >(pstart,pstop),*this);
502 else
503 Kokkos::parallel_reduce(Kokkos::RangePolicy<DeviceType, TagUpdateMove<DIM,SURF,-1> >(pstart,pstop),*this,reduce);
504 copymode = 0;
505
506 Kokkos::deep_copy(h_scalars,d_scalars);
507
508 particle_kk->modify(Device,PARTICLE_MASK);
509 d_particles = t_particle_1d(); // destroy reference to reduce memory use
510
511 k_mlist.modify_device();
512
513 // END of pstart/pstop loop advecting all particles
514
515 nmigrate = h_nmigrate();
516
517 int error_flag;
518
519 if (!sparta->kokkos->need_atomics || sparta->kokkos->atomic_reduction) {
520 ntouch_one = h_ntouch_one();
521 nexit_one = h_nexit_one();
522 nboundary_one = h_nboundary_one();
523 ncomm_one = h_ncomm_one();
524 nscheck_one = h_nscheck_one();
525 nscollide_one = h_nscollide_one();
526 surf->nreact_one = h_nreact_one();
527 nstuck = h_nstuck();
528 naxibad = h_naxibad();
529 } else {
530 ntouch_one += reduce.ntouch_one ;
531 nexit_one += reduce.nexit_one ;
532 nboundary_one += reduce.nboundary_one;
533 ncomm_one += reduce.ncomm_one ;
534 nscheck_one += reduce.nscheck_one ;
535 nscollide_one += reduce.nscollide_one;
536 surf->nreact_one += reduce.nreact_one ;
537 nstuck += reduce.nstuck ;
538 naxibad += reduce.nstuck ;
539 }
540
541 entryexit = h_entryexit();
542
543 error_flag = h_error_flag();
544
545 if (error_flag) {
546 char str[128];
547 sprintf(str,
548 "Particle being sent to self proc "
549 "on step " BIGINT_FORMAT,
550 update->ntimestep);
551 error->one(FLERR,str);
552 }
553
554 // if gridcut >= 0.0, check if another iteration of move is required
555 // only the case if some particle flag = PENTRY/PEXIT
556 // in which case perform particle migration
557 // if not, move is done and final particle comm will occur in run()
558 // if iterating, reset pstart/pstop and extend migration list if necessary
559
560 if (grid->cutoff < 0.0) break;
561
562 timer->stamp(TIME_MOVE);
563 MPI_Allreduce(&entryexit,&any_entryexit,1,MPI_INT,MPI_MAX,world);
564 timer->stamp();
565 if (any_entryexit) {
566 if (nmigrate) {
567 k_mlist_small = Kokkos::subview(k_mlist,std::make_pair(0,nmigrate));
568 k_mlist_small.sync_host();
569 }
570 auto mlist_small = k_mlist_small.h_view.data();
571 timer->stamp(TIME_MOVE);
572 pstart = ((CommKokkos*)comm)->migrate_particles(nmigrate,mlist_small,k_mlist_small.d_view);
573 timer->stamp(TIME_COMM);
574 pstop = particle->nlocal;
575 if (pstop-pstart > maxmigrate) {
576 maxmigrate = pstop-pstart;
577 memoryKK->destroy_kokkos(k_mlist,mlist);
578 memoryKK->create_kokkos(k_mlist,mlist,maxmigrate,"particle:mlist");
579 }
580 } else break;
581
582 // END of single move/migrate iteration
583
584 }
585
586 // END of all move/migrate iterations
587
588 particle->sorted = 0;
589 particle_kk->sorted_kk = 0;
590
591 // accumulate running totals
592
593 niterate_running += niterate;
594 nmove_running += nlocal;
595 ntouch_running += ntouch_one;
596 ncomm_running += ncomm_one;
597 nboundary_running += nboundary_one;
598 nexit_running += nexit_one;
599 nscheck_running += nscheck_one;
600 nscollide_running += nscollide_one;
601 surf->nreact_running += surf->nreact_one;
602
603 if (nsurf_tally) {
604 for (int m = 0; m < nsurf_tally; m++) {
605 ComputeSurfKokkos* compute_surf_kk = (ComputeSurfKokkos*)(slist_active[m]);
606 compute_surf_kk->post_surf_tally();
607 }
608 }
609
610 if (nboundary_tally) {
611 for (int m = 0; m < nboundary_tally; m++) {
612 ComputeBoundaryKokkos* compute_boundary_kk = (ComputeBoundaryKokkos*)(blist_active[m]);
613 compute_boundary_kk->post_boundary_tally();
614 }
615 }
616 }
617
618 /* ---------------------------------------------------------------------- */
619
620 template<int DIM, int SURF, int ATOMIC_REDUCTION>
621 KOKKOS_INLINE_FUNCTION
operator ()(TagUpdateMove<DIM,SURF,ATOMIC_REDUCTION>,const int & i) const622 void UpdateKokkos::operator()(TagUpdateMove<DIM,SURF,ATOMIC_REDUCTION>, const int &i) const {
623 UPDATE_REDUCE reduce;
624 this->template operator()<DIM,SURF,ATOMIC_REDUCTION>(TagUpdateMove<DIM,SURF,ATOMIC_REDUCTION>(), i, reduce);
625 }
626
627 /*-----------------------------------------------------------------------------*/
628
629 template<int DIM, int SURF, int ATOMIC_REDUCTION>
630 KOKKOS_INLINE_FUNCTION
operator ()(TagUpdateMove<DIM,SURF,ATOMIC_REDUCTION>,const int & i,UPDATE_REDUCE & reduce) const631 void UpdateKokkos::operator()(TagUpdateMove<DIM,SURF,ATOMIC_REDUCTION>, const int &i, UPDATE_REDUCE &reduce) const {
632 if (d_error_flag()) return;
633
634 // int m;
635 bool hitflag;
636 int icell,icell_original,outface,bflag,nflag,pflag,itmp;
637 int side,minsurf,nsurf,cflag,isurf,exclude,stuck_iterate;
638 double dtremain,frac,newfrac,param,minparam,rnew,dtsurf,tc,tmp;
639 double xnew[3],xhold[3],xc[3],vc[3],minxc[3],minvc[3];
640 double *x,*v;
641 Surf::Tri *tri;
642 Surf::Line *line;
643 int reaction; // not yet used
644
645 Particle::OnePart &particle_i = d_particles[i];
646 pflag = particle_i.flag;
647
648 Particle::OnePart iorig;
649 Particle::OnePart *ipart,*jpart;
650 jpart = NULL;
651
652 // received from another proc and move is done
653 // if first iteration, PDONE is from a previous step,
654 // set pflag to PKEEP so move the particle on this step
655 // else do nothing
656
657 if (pflag == PDONE) {
658 pflag = particle_i.flag = PKEEP;
659 if (niterate > 1) return;
660 }
661
662 x = particle_i.x;
663 v = particle_i.v;
664 exclude = -1;
665
666 // for 2d and axisymmetry only
667 // xnew,xc passed to geometry routines which use or set z component
668
669 if (DIM < 3) xnew[2] = xc[2] = 0.0;
670
671 // apply moveperturb() to PKEEP and PINSERT since are computing xnew
672 // not to PENTRY,PEXIT since are just re-computing xnew of sender
673 // set xnew[2] to linear move for axisymmetry, will be remapped later
674 // let pflag = PEXIT persist to check during axisymmetric cell crossing
675
676 if (DIM < 3) xnew[2] = 0.0;
677 if (pflag == PKEEP) {
678 dtremain = dt;
679 xnew[0] = x[0] + dtremain*v[0];
680 xnew[1] = x[1] + dtremain*v[1];
681 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
682 if (gravity_3d_flag) gravity3d(dtremain,xnew,v);
683 else if (gravity_2d_flag) gravity2d(dtremain,xnew,v);
684 } else if (pflag == PINSERT) {
685 dtremain = particle_i.dtremain;
686 xnew[0] = x[0] + dtremain*v[0];
687 xnew[1] = x[1] + dtremain*v[1];
688 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
689 if (gravity_3d_flag) gravity3d(dtremain,xnew,v);
690 else if (gravity_2d_flag) gravity2d(dtremain,xnew,v);
691 } else if (pflag == PENTRY) {
692 icell = particle_i.icell;
693 if (d_cells[icell].nsplit > 1) {
694 if (DIM == 3 && SURF) icell = split3d(icell,x);
695 if (DIM < 3 && SURF) icell = split2d(icell,x);
696 particle_i.icell = icell;
697 }
698 dtremain = particle_i.dtremain;
699 xnew[0] = x[0] + dtremain*v[0];
700 xnew[1] = x[1] + dtremain*v[1];
701 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
702 } else if (pflag == PEXIT) {
703 dtremain = particle_i.dtremain;
704 xnew[0] = x[0] + dtremain*v[0];
705 xnew[1] = x[1] + dtremain*v[1];
706 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
707 } else if (pflag >= PSURF) {
708 dtremain = particle_i.dtremain;
709 xnew[0] = x[0] + dtremain*v[0];
710 xnew[1] = x[1] + dtremain*v[1];
711 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
712 if (pflag > PSURF) exclude = pflag - PSURF - 1;
713 }
714
715 particle_i.flag = PKEEP;
716 icell = particle_i.icell;
717 double* lo = d_cells[icell].lo;
718 double* hi = d_cells[icell].hi;
719 cellint* neigh = d_cells[icell].neigh;
720 int nmask = d_cells[icell].nmask;
721 stuck_iterate = 0;
722 if (ATOMIC_REDUCTION == 1)
723 Kokkos::atomic_fetch_add(&d_ntouch_one(),1);
724 else if (ATOMIC_REDUCTION == 0)
725 d_ntouch_one()++;
726 else
727 reduce.ntouch_one++;
728
729 // advect one particle from cell to cell and thru surf collides til done
730
731 while (1) {
732
733 #ifdef MOVE_DEBUG
734 if (DIM == 3) {
735 if (ntimestep == MOVE_DEBUG_STEP &&
736 (MOVE_DEBUG_ID == d_particles[i].id ||
737 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
738 printf("PARTICLE %d %ld: %d %d: %d: x %g %g %g: xnew %g %g %g: %d "
739 CELLINT_FORMAT ": lo %g %g %g: hi %g %g %g: DTR %g\n",
740 me,ntimestep,i,d_particles[i].id,
741 d_cells[icell].nsurf,
742 x[0],x[1],x[2],xnew[0],xnew[1],xnew[2],
743 icell,d_cells[icell].id,
744 lo[0],lo[1],lo[2],hi[0],hi[1],hi[2],dtremain);
745 }
746 if (DIM == 2) {
747 if (ntimestep == MOVE_DEBUG_STEP &&
748 (MOVE_DEBUG_ID == d_particles[i].id ||
749 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
750 printf("PARTICLE %d %ld: %d %d: %d: x %g %g: xnew %g %g: %d "
751 CELLINT_FORMAT ": lo %g %g: hi %g %g: DTR: %g\n",
752 me,ntimestep,i,d_particles[i].id,
753 d_cells[icell].nsurf,
754 x[0],x[1],xnew[0],xnew[1],
755 icell,d_cells[icell].id,
756 lo[0],lo[1],hi[0],hi[1],dtremain);
757 }
758 if (DIM == 1) {
759 if (ntimestep == MOVE_DEBUG_STEP &&
760 (MOVE_DEBUG_ID == d_particles[i].id ||
761 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
762 printf("PARTICLE %d %ld: %d %d: %d: x %g %g: xnew %g %g: %d "
763 CELLINT_FORMAT ": lo %g %g: hi %g %g: DTR: %g\n",
764 me,ntimestep,i,d_particles[i].id,
765 d_cells[icell].nsurf,
766 x[0],x[1],xnew[0],sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]),
767 icell,d_cells[icell].id,
768 lo[0],lo[1],hi[0],hi[1],dtremain);
769 }
770 #endif
771
772 // check if particle crosses any cell face
773 // frac = fraction of move completed before hitting cell face
774 // this section should be as efficient as possible,
775 // since most particles won't do anything else
776 // axisymmetric cell face crossings:
777 // use linear xnew to check vertical faces
778 // must always check move against curved lower y face of cell
779 // use remapped rnew to check horizontal lines
780 // for y faces, if pflag = PEXIT, particle was just received
781 // from another proc and is exiting this cell from face:
782 // axi_horizontal_line() will not detect correct crossing,
783 // so set frac and outface directly to move into adjacent cell,
784 // then unset pflag so not checked again for this particle
785
786 outface = INTERIOR;
787 frac = 1.0;
788
789 if (xnew[0] < lo[0]) {
790 frac = (lo[0]-x[0]) / (xnew[0]-x[0]);
791 outface = XLO;
792 } else if (xnew[0] >= hi[0]) {
793 frac = (hi[0]-x[0]) / (xnew[0]-x[0]);
794 outface = XHI;
795 }
796
797 if (DIM != 1) {
798 if (xnew[1] < lo[1]) {
799 newfrac = (lo[1]-x[1]) / (xnew[1]-x[1]);
800 if (newfrac < frac) {
801 frac = newfrac;
802 outface = YLO;
803 }
804 } else if (xnew[1] >= hi[1]) {
805 newfrac = (hi[1]-x[1]) / (xnew[1]-x[1]);
806 if (newfrac < frac) {
807 frac = newfrac;
808 outface = YHI;
809 }
810 }
811 }
812
813 if (DIM == 1) {
814 if (x[1] == lo[1] && (pflag == PEXIT || v[1] < 0.0)) {
815 frac = 0.0;
816 outface = YLO;
817 } else if (GeometryKokkos::
818 axi_horizontal_line(dtremain,x,v,lo[1],itmp,tc,tmp)) {
819 newfrac = tc/dtremain;
820 if (newfrac < frac) {
821 frac = newfrac;
822 outface = YLO;
823 }
824 }
825
826 if (x[1] == hi[1] && (pflag == PEXIT || v[1] > 0.0)) {
827 frac = 0.0;
828 outface = YHI;
829 } else {
830 rnew = sqrt(xnew[1]*xnew[1] + xnew[2]*xnew[2]);
831 if (rnew >= hi[1]) {
832 if (GeometryKokkos::
833 axi_horizontal_line(dtremain,x,v,hi[1],itmp,tc,tmp)) {
834 newfrac = tc/dtremain;
835 if (newfrac < frac) {
836 frac = newfrac;
837 outface = YHI;
838 }
839 }
840 }
841 }
842
843 pflag = 0;
844 }
845
846 if (DIM == 3) {
847 if (xnew[2] < lo[2]) {
848 newfrac = (lo[2]-x[2]) / (xnew[2]-x[2]);
849 if (newfrac < frac) {
850 frac = newfrac;
851 outface = ZLO;
852 }
853 } else if (xnew[2] >= hi[2]) {
854 newfrac = (hi[2]-x[2]) / (xnew[2]-x[2]);
855 if (newfrac < frac) {
856 frac = newfrac;
857 outface = ZHI;
858 }
859 }
860 }
861
862 #ifdef MOVE_DEBUG
863 if (ntimestep == MOVE_DEBUG_STEP &&
864 (MOVE_DEBUG_ID == d_particles[i].id ||
865 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX))) {
866 if (outface != INTERIOR)
867 printf(" OUTFACE %d out: %d %d, frac %g\n",
868 outface,grid_kk_copy.obj.neigh_decode(nmask,outface),
869 neigh[outface],frac);
870 else
871 printf(" INTERIOR %d %d\n",outface,INTERIOR);
872 }
873 #endif
874
875 // START of code specific to surfaces
876
877 if (SURF) {
878
879 // skip surf checks if particle flagged as EXITing this cell
880 // then unset pflag so not checked again for this particle
881
882 nsurf = d_cells[icell].nsurf;
883 if (pflag == PEXIT) {
884 nsurf = 0;
885 pflag = 0;
886 }
887
888 if (ATOMIC_REDUCTION == 1)
889 Kokkos::atomic_fetch_add(&d_nscheck_one(),nsurf);
890 else if (ATOMIC_REDUCTION == 0)
891 d_nscheck_one() += nsurf;
892 else
893 reduce.nscheck_one += nsurf;
894
895 if (nsurf) {
896
897 // particle crosses cell face, reset xnew exactly on face of cell
898 // so surface check occurs only for particle path within grid cell
899 // xhold = saved xnew so can restore below if no surf collision
900
901 if (outface != INTERIOR) {
902 xhold[0] = xnew[0];
903 xhold[1] = xnew[1];
904 if (DIM != 2) xhold[2] = xnew[2];
905
906 xnew[0] = x[0] + frac*(xnew[0]-x[0]);
907 xnew[1] = x[1] + frac*(xnew[1]-x[1]);
908 if (DIM != 2) xnew[2] = x[2] + frac*(xnew[2]-x[2]);
909
910 if (outface == XLO) xnew[0] = lo[0];
911 else if (outface == XHI) xnew[0] = hi[0];
912 else if (outface == YLO) xnew[1] = lo[1];
913 else if (outface == YHI) xnew[1] = hi[1];
914 else if (outface == ZLO) xnew[2] = lo[2];
915 else if (outface == ZHI) xnew[2] = hi[2];
916 }
917
918 // for axisymmetric, dtsurf = time that particle stays in cell
919 // used as arg to axi_line_intersect()
920
921 if (DIM == 1) {
922 if (outface == INTERIOR) dtsurf = dtremain;
923 else dtsurf = dtremain * frac;
924 }
925
926 // check for collisions with triangles or lines in cell
927 // find 1st surface hit via minparam
928 // skip collisions with previous surf, but not for axisymmetric
929 // not considered collision if 2 params are tied and one INSIDE surf
930 // if collision occurs, perform collision with surface model
931 // reset x,v,xnew,dtremain and continue single particle trajectory
932
933 cflag = 0;
934 minparam = 2.0;
935 auto csurfs_begin = d_csurfs.row_map(icell);
936
937 for (int m = 0; m < nsurf; m++) {
938 isurf = d_csurfs.entries(csurfs_begin + m);
939
940 if (DIM > 1) {
941 if (isurf == exclude) continue;
942 }
943 if (DIM == 3) {
944 tri = &d_tris[isurf];
945 hitflag = GeometryKokkos::
946 line_tri_intersect(x,xnew,
947 tri->p1,tri->p2,
948 tri->p3,tri->norm,xc,param,side);
949 }
950 if (DIM == 2) {
951 line = &d_lines[isurf];
952 hitflag = GeometryKokkos::
953 line_line_intersect(x,xnew,
954 line->p1,line->p2,
955 line->norm,xc,param,side);
956 }
957 if (DIM == 1) {
958 line = &d_lines[isurf];
959 hitflag = GeometryKokkos::
960 axi_line_intersect(dtsurf,x,v,outface,lo,hi,
961 line->p1,line->p2,
962 line->norm,exclude == isurf,
963 xc,vc,param,side);
964 }
965
966 #ifdef MOVE_DEBUG
967 if (DIM == 3) {
968 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
969 (MOVE_DEBUG_ID == d_particles[i].id ||
970 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
971 printf("SURF COLLIDE: %d %d %d %d: "
972 "P1 %g %g %g: P2 %g %g %g: "
973 "T1 %g %g %g: T2 %g %g %g: T3 %g %g %g: "
974 "TN %g %g %g: XC %g %g %g: "
975 "Param %g: Side %d\n",
976 MOVE_DEBUG_INDEX,icell,nsurf,isurf,
977 x[0],x[1],x[2],xnew[0],xnew[1],xnew[2],
978 tri->p1[0],tri->p1[1],tri->p1[2],
979 tri->p2[0],tri->p2[1],tri->p2[2],
980 tri->p3[0],tri->p3[1],tri->p3[2],
981 tri->norm[0],tri->norm[1],tri->norm[2],
982 xc[0],xc[1],xc[2],param,side);
983 }
984 if (DIM == 2) {
985 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
986 (MOVE_DEBUG_ID == d_particles[i].id ||
987 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
988 printf("SURF COLLIDE: %d %d %d %d: P1 %g %g: P2 %g %g: "
989 "L1 %g %g: L2 %g %g: LN %g %g: XC %g %g: "
990 "Param %g: Side %d\n",
991 MOVE_DEBUG_INDEX,icell,nsurf,isurf,
992 x[0],x[1],xnew[0],xnew[1],
993 line->p1[0],line->p1[1],line->p2[0],line->p2[1],
994 line->norm[0],line->norm[1],
995 xc[0],xc[1],param,side);
996 }
997 if (DIM == 1) {
998 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
999 (MOVE_DEBUG_ID == d_particles[i].id ||
1000 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
1001 printf("SURF COLLIDE %d %ld: %d %d %d %d: P1 %g %g: P2 %g %g: "
1002 "L1 %g %g: L2 %g %g: LN %g %g: XC %g %g: "
1003 "VC %g %g %g: Param %g: Side %d\n",
1004 hitflag,ntimestep,MOVE_DEBUG_INDEX,icell,nsurf,isurf,
1005 x[0],x[1],
1006 xnew[0],sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]),
1007 line->p1[0],line->p1[1],line->p2[0],line->p2[1],
1008 line->norm[0],line->norm[1],
1009 xc[0],xc[1],vc[0],vc[1],vc[2],param,side);
1010 double edge1[3],edge2[3],xfinal[3],cross[3];
1011 MathExtra::sub3(line->p2,line->p1,edge1);
1012 MathExtra::sub3(x,line->p1,edge2);
1013 MathExtra::cross3(edge2,edge1,cross);
1014 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
1015 MOVE_DEBUG_ID == d_particles[i].id)
1016 printf("CROSSSTART %g %g %g\n",cross[0],cross[1],cross[2]);
1017 xfinal[0] = xnew[0];
1018 xfinal[1] = sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]);
1019 xfinal[2] = 0.0;
1020 MathExtra::sub3(xfinal,line->p1,edge2);
1021 MathExtra::cross3(edge2,edge1,cross);
1022 if (hitflag && ntimestep == MOVE_DEBUG_STEP &&
1023 MOVE_DEBUG_ID == d_particles[i].id)
1024 printf("CROSSFINAL %g %g %g\n",cross[0],cross[1],cross[2]);
1025 }
1026 #endif
1027
1028 if (hitflag && param < minparam && side == OUTSIDE) {
1029
1030 // NOTE: these were the old checks
1031 // think it is now sufficient to test for particle
1032 // in an INSIDE cell in fix grid/check
1033
1034 //if (hitflag && side != ONSURF2OUT && param <= minparam)
1035
1036 // this if test is to avoid case where particle
1037 // previously hit 1 of 2 (or more) touching angled surfs at
1038 // common edge/corner, on this iteration first surf
1039 // is excluded, but others may be hit on inside:
1040 // param will be epsilon and exclude must be set
1041 // skip the hits of other touching surfs
1042
1043 //if (side == INSIDE && param < EPSPARAM && exclude >= 0)
1044 // continue;
1045
1046 // this if test is to avoid case where particle
1047 // hits 2 touching angled surfs at common edge/corner
1048 // from far away:
1049 // param is same, but hits one on outside, one on inside
1050 // only keep surf hit on outside
1051
1052 //if (param == minparam && side == INSIDE) continue;
1053
1054 cflag = 1;
1055 minparam = param;
1056 // minside = side;
1057 minsurf = isurf;
1058 minxc[0] = xc[0];
1059 minxc[1] = xc[1];
1060 if (DIM == 3) minxc[2] = xc[2];
1061 if (DIM == 1) {
1062 minvc[1] = vc[1];
1063 minvc[2] = vc[2];
1064 }
1065 }
1066
1067 } // END of for loop over surfs
1068
1069 // tri/line = surf that particle hit first
1070
1071 if (cflag) {
1072 if (DIM == 3) tri = &d_tris[minsurf];
1073 if (DIM != 3) line = &d_lines[minsurf];
1074
1075 // set x to collision point
1076 // if axisymmetric, set v to remapped velocity at collision pt
1077
1078 x[0] = minxc[0];
1079 x[1] = minxc[1];
1080 if (DIM == 3) x[2] = minxc[2];
1081 if (DIM == 1) {
1082 v[1] = minvc[1];
1083 v[2] = minvc[2];
1084 }
1085
1086 // perform surface collision using surface collision model
1087 // surface chemistry may destroy particle or create new one
1088 // must update particle's icell to current icell so that
1089 // if jpart is created, it will be added to correct cell
1090 // if jpart, add new particle to this iteration via pstop++
1091 // tally surface statistics if requested using iorig
1092
1093 ipart = &particle_i;
1094 ipart->icell = icell;
1095 dtremain *= 1.0 - minparam*frac;
1096
1097 if (nsurf_tally)
1098 iorig = particle_i;
1099 int n = DIM == 3 ? tri->isc : line->isc;
1100 int sc_type = sc_type_list[n];
1101 int m = sc_map[n];
1102
1103 if (DIM == 3)
1104 if (sc_type == 0)
1105 jpart = sc_kk_specular_copy[m].obj.
1106 collide_kokkos(ipart,tri->norm,dtremain,tri->isr,reaction);/////
1107 else if (sc_type == 1)
1108 jpart = sc_kk_diffuse_copy[m].obj.
1109 collide_kokkos(ipart,tri->norm,dtremain,tri->isr,reaction);/////
1110 else if (sc_type == 2)
1111 jpart = sc_kk_vanish_copy[m].obj.
1112 collide_kokkos(ipart,tri->norm,dtremain,tri->isr,reaction);/////
1113 else if (sc_type == 3)
1114 jpart = sc_kk_piston_copy[m].obj.
1115 collide_kokkos(ipart,tri->norm,dtremain,tri->isr,reaction);/////
1116 else if (sc_type == 4)
1117 jpart = sc_kk_transparent_copy[m].obj.
1118 collide_kokkos(ipart,tri->norm,dtremain,tri->isr,reaction);/////
1119 if (DIM != 3)
1120 if (sc_type == 0)
1121 jpart = sc_kk_specular_copy[m].obj.
1122 collide_kokkos(ipart,line->norm,dtremain,line->isr,reaction);////
1123 else if (sc_type == 1)
1124 jpart = sc_kk_diffuse_copy[m].obj.
1125 collide_kokkos(ipart,line->norm,dtremain,line->isr,reaction);////
1126 else if (sc_type == 2)
1127 jpart = sc_kk_vanish_copy[m].obj.
1128 collide_kokkos(ipart,line->norm,dtremain,line->isr,reaction);////
1129 else if (sc_type == 3)
1130 jpart = sc_kk_piston_copy[m].obj.
1131 collide_kokkos(ipart,line->norm,dtremain,line->isr,reaction);////
1132 else if (sc_type == 4)
1133 jpart = sc_kk_transparent_copy[m].obj.
1134 collide_kokkos(ipart,line->norm,dtremain,line->isr,reaction);////
1135
1136 ////Need to error out for now if surface reactions create (or destroy?) particles////
1137 //if (jpart) {
1138 // particles = particle->particles;
1139 // x = particle_i.x;
1140 // v = particle_i.v;
1141 // jpart->flag = PSURF + 1 + minsurf;
1142 // jpart->dtremain = dtremain;
1143 // jpart->weight = particle_i.weight;
1144 // pstop++;
1145 //}
1146
1147 if (nsurf_tally)
1148 for (m = 0; m < nsurf_tally; m++)
1149 slist_active_copy[m].obj.
1150 surf_tally_kk<ATOMIC_REDUCTION>(minsurf,icell,reaction,&iorig,ipart,jpart);
1151
1152 // stuck_iterate = consecutive iterations particle is immobile
1153
1154 if (minparam == 0.0) stuck_iterate++;
1155 else stuck_iterate = 0;
1156
1157 // reset post-bounce xnew
1158
1159 xnew[0] = x[0] + dtremain*v[0];
1160 xnew[1] = x[1] + dtremain*v[1];
1161 if (DIM != 2) xnew[2] = x[2] + dtremain*v[2];
1162
1163 exclude = minsurf;
1164 if (ATOMIC_REDUCTION == 1)
1165 Kokkos::atomic_fetch_add(&d_nscollide_one(),1);
1166 else if (ATOMIC_REDUCTION == 0)
1167 d_nscollide_one()++;
1168 else
1169 reduce.nscollide_one++;
1170
1171 #ifdef MOVE_DEBUG
1172 if (DIM == 3) {
1173 if (ntimestep == MOVE_DEBUG_STEP &&
1174 (MOVE_DEBUG_ID == d_particles[i].id ||
1175 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
1176 printf("POST COLLISION %d: %g %g %g: %g %g %g: %g %g %g\n",
1177 MOVE_DEBUG_INDEX,
1178 x[0],x[1],x[2],xnew[0],xnew[1],xnew[2],
1179 minparam,frac,dtremain);
1180 }
1181 if (DIM == 2) {
1182 if (ntimestep == MOVE_DEBUG_STEP &&
1183 (MOVE_DEBUG_ID == d_particles[i].id ||
1184 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
1185 printf("POST COLLISION %d: %g %g: %g %g: %g %g %g\n",
1186 MOVE_DEBUG_INDEX,
1187 x[0],x[1],xnew[0],xnew[1],
1188 minparam,frac,dtremain);
1189 }
1190 if (DIM == 1) {
1191 if (ntimestep == MOVE_DEBUG_STEP &&
1192 (MOVE_DEBUG_ID == d_particles[i].id ||
1193 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
1194 printf("POST COLLISION %d: %g %g: %g %g: vel %g %g %g: %g %g %g\n",
1195 MOVE_DEBUG_INDEX,
1196 x[0],x[1],
1197 xnew[0],sqrt(xnew[1]*xnew[1]+xnew[2]*xnew[2]),
1198 v[0],v[1],v[2],
1199 minparam,frac,dtremain);
1200 }
1201 #endif
1202
1203 // if ipart = NULL, particle discarded due to surface chem
1204 // else if particle not stuck, continue advection while loop
1205 // if stuck, mark for DISCARD, and drop out of SURF code
1206
1207 if (ipart == NULL) particle_i.flag = PDISCARD;
1208 else if (stuck_iterate < MAXSTUCK) continue;
1209 else {
1210 particle_i.flag = PDISCARD;
1211 if (ATOMIC_REDUCTION == 1)
1212 Kokkos::atomic_fetch_add(&d_nstuck(),1);
1213 else if (ATOMIC_REDUCTION == 0)
1214 d_nstuck()++;
1215 else
1216 reduce.nstuck++;
1217 }
1218
1219 } // END of cflag if section that performed collision
1220
1221 // no collision, so restore saved xnew if changed it above
1222
1223 if (outface != INTERIOR) {
1224 xnew[0] = xhold[0];
1225 xnew[1] = xhold[1];
1226 if (DIM != 2) xnew[2] = xhold[2];
1227 }
1228
1229 } // END of if test for any surfs in this cell
1230 } // END of code specific to surfaces
1231
1232 // break from advection loop if discarding particle
1233
1234 if (particle_i.flag == PDISCARD) break;
1235
1236
1237 // no cell crossing
1238 // set final particle position to xnew, then break from advection loop
1239 // for axisymmetry, must first remap linear xnew and v
1240 // for axisymmetry, check if final particle position is within cell
1241 // can be rare epsilon round-off cases where particle ends up outside
1242 // of final cell curved surf when move logic thinks it is inside
1243 // example is when Geom::axi_horizontal_line() says no crossing of cell edge
1244 // but axi_remap() puts particle outside the cell
1245 // in this case, just DISCARD particle and tally it to naxibad
1246 // if migrating to another proc,
1247 // flag as PDONE so new proc won't move it more on this step
1248
1249 if (outface == INTERIOR) {
1250 if (DIM == 1) axi_remap(xnew,v);
1251 x[0] = xnew[0];
1252 x[1] = xnew[1];
1253 if (DIM == 3) x[2] = xnew[2];
1254 if (DIM == 1) {
1255 if (x[1] < lo[1] || x[1] > hi[1]) {
1256 particle_i.flag = PDISCARD;
1257 if (ATOMIC_REDUCTION == 1)
1258 Kokkos::atomic_fetch_add(&d_naxibad(),1);
1259 else if (ATOMIC_REDUCTION == 0)
1260 d_naxibad()++;
1261 else
1262 reduce.naxibad++;
1263 break;
1264 }
1265 }
1266 if (d_cells[icell].proc != me) particle_i.flag = PDONE;
1267 break;
1268 }
1269
1270 // particle crosses cell face
1271 // decrement dtremain in case particle is passed to another proc
1272 // for axisymmetry, must then remap linear x and v
1273 // reset particle x to be exactly on cell face
1274 // for axisymmetry, must reset xnew for next iteration since v changed
1275
1276 dtremain *= 1.0-frac;
1277 exclude = -1;
1278
1279 x[0] += frac * (xnew[0]-x[0]);
1280 x[1] += frac * (xnew[1]-x[1]);
1281 if (DIM != 2) x[2] += frac * (xnew[2]-x[2]);
1282 if (DIM == 1) axi_remap(x,v);
1283
1284 if (outface == XLO) x[0] = lo[0];
1285 else if (outface == XHI) x[0] = hi[0];
1286 else if (outface == YLO) x[1] = lo[1];
1287 else if (outface == YHI) x[1] = hi[1];
1288 else if (outface == ZLO) x[2] = lo[2];
1289 else if (outface == ZHI) x[2] = hi[2];
1290
1291 if (DIM == 1) {
1292 xnew[0] = x[0] + dtremain*v[0];
1293 xnew[1] = x[1] + dtremain*v[1];
1294 xnew[2] = x[2] + dtremain*v[2];
1295 }
1296
1297 // nflag = type of neighbor cell: child, parent, unknown, boundary
1298 // if parent, use id_find_child to identify child cell
1299 // result can be -1 for unknown cell, occurs when:
1300 // (a) particle hits face of ghost child cell
1301 // (b) the ghost cell extends beyond ghost halo
1302 // (c) cell on other side of face is a parent
1303 // (d) its child, which the particle is in, is entirely beyond my halo
1304 // if new cell is child and surfs exist, check if a split cell
1305
1306 nflag = grid_kk_copy.obj.neigh_decode(nmask,outface);
1307 icell_original = icell;
1308
1309 if (nflag == NCHILD) {
1310 icell = neigh[outface];
1311 if (DIM == 3 && SURF) {
1312 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1313 icell = split3d(icell,x);
1314 }
1315 if (DIM < 3 && SURF) {
1316 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1317 icell = split2d(icell,x);
1318 }
1319 } else if (nflag == NPARENT) {
1320 auto pcell = &d_pcells[neigh[outface]];
1321 icell = grid_kk_copy.obj.id_find_child(pcell->id,d_cells[icell].level,
1322 pcell->lo,pcell->hi,x);
1323 if (icell >= 0) {
1324 if (DIM == 3 && SURF) {
1325 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1326 icell = split3d(icell,x);
1327 }
1328 if (DIM < 3 && SURF) {
1329 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1330 icell = split2d(icell,x);
1331 }
1332 }
1333 } else if (nflag == NUNKNOWN) icell = -1;
1334
1335 // neighbor cell is global boundary
1336 // tally boundary stats if requested using iorig
1337 // collide() updates x,v,xnew as needed due to boundary interaction
1338 // may also update dtremain (piston BC)
1339 // for axisymmetric, must recalculate xnew since v may have changed
1340 // surface chemistry may destroy particle or create new one
1341 // if jpart, add new particle to this iteration via pstop++
1342 // OUTFLOW: exit with particle flag = PDISCARD
1343 // PERIODIC: new cell via same logic as above for child/parent/unknown
1344 // OTHER = reflected particle stays in same grid cell
1345
1346 else {
1347 ipart = &particle_i;
1348
1349 Particle::OnePart iorig;
1350 if (nboundary_tally)
1351 memcpy(&iorig,&particle_i,sizeof(Particle::OnePart));
1352
1353 Particle::OnePart* ipart = &particle_i;
1354 lo = d_cells[icell].lo;
1355 hi = d_cells[icell].hi;
1356 if (domain_kk_copy.obj.bflag[outface] == SURFACE) {
1357 // treat global boundary as a surface
1358 // particle velocity is changed by surface collision model
1359 // dtremain may be changed by collision model
1360 // reset all components of xnew, in case dtremain changed
1361 // if axisymmetric, caller will reset again, including xnew[2]
1362
1363 int n = domain_kk_copy.obj.surf_collide[outface];
1364 int sc_type = sc_type_list[n];
1365 int m = sc_map[n];
1366
1367 if (sc_type == 0)
1368 jpart = sc_kk_specular_copy[m].obj.
1369 collide_kokkos(ipart,domain_kk_copy.obj.norm[outface],dtremain,domain_kk_copy.obj.surf_react[outface],reaction);/////
1370 else if (sc_type == 1)
1371 jpart = sc_kk_diffuse_copy[m].obj.
1372 collide_kokkos(ipart,domain_kk_copy.obj.norm[outface],dtremain,domain_kk_copy.obj.surf_react[outface],reaction);/////
1373 else if (sc_type == 2)
1374 jpart = sc_kk_vanish_copy[m].obj.
1375 collide_kokkos(ipart,domain_kk_copy.obj.norm[outface],dtremain,domain_kk_copy.obj.surf_react[outface],reaction);/////
1376 else if (sc_type == 3)
1377 jpart = sc_kk_piston_copy[m].obj.
1378 collide_kokkos(ipart,domain_kk_copy.obj.norm[outface],dtremain,domain_kk_copy.obj.surf_react[outface],reaction);/////
1379 else if (sc_type == 4)
1380 jpart = sc_kk_transparent_copy[m].obj.
1381 collide_kokkos(ipart,domain_kk_copy.obj.norm[outface],dtremain,domain_kk_copy.obj.surf_react[outface],reaction);/////
1382
1383 if (ipart) {
1384 double *x = ipart->x;
1385 double *v = ipart->v;
1386 xnew[0] = x[0] + dtremain*v[0];
1387 xnew[1] = x[1] + dtremain*v[1];
1388 if (domain_kk_copy.obj.dimension == 3) xnew[2] = x[2] + dtremain*v[2];
1389 }
1390 bflag = SURFACE;
1391 } else {
1392 bflag = domain_kk_copy.obj.collide_kokkos(ipart,outface,lo,hi,xnew/*,dtremain*/,reaction);
1393 }
1394
1395 //if (jpart) {
1396 // particles = particle->particles;
1397 // x = particle_i.x;
1398 // v = particle_i.v;
1399 //}
1400
1401 if (nboundary_tally)
1402 for (int m = 0; m < nboundary_tally; m++)
1403 blist_active_copy[m].obj.
1404 boundary_tally_kk<ATOMIC_REDUCTION>(outface,bflag,reaction,&iorig,ipart,jpart,domain_kk_copy.obj.norm[outface]);
1405
1406 if (DIM == 1) {
1407 xnew[0] = x[0] + dtremain*v[0];
1408 xnew[1] = x[1] + dtremain*v[1];
1409 xnew[2] = x[2] + dtremain*v[2];
1410 }
1411
1412 if (bflag == OUTFLOW) {
1413 particle_i.flag = PDISCARD;
1414 if (ATOMIC_REDUCTION == 1)
1415 Kokkos::atomic_fetch_add(&d_nexit_one(),1);
1416 else if (ATOMIC_REDUCTION == 0)
1417 d_nexit_one()++;
1418 else
1419 reduce.nexit_one++;
1420 break;
1421 } else if (bflag == PERIODIC) {
1422 if (nflag == NPBCHILD) {
1423 icell = neigh[outface];
1424 if (DIM == 3 && SURF) {
1425 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1426 icell = split3d(icell,x);
1427 }
1428 if (DIM < 3 && SURF) {
1429 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1430 icell = split2d(icell,x);
1431 }
1432 } else if (nflag == NPBPARENT) {
1433 auto pcell = &d_pcells[neigh[outface]];
1434 icell = grid_kk_copy.obj.id_find_child(pcell->id,d_cells[icell].level,
1435 pcell->lo,pcell->hi,x);
1436 if (icell >= 0) {
1437 if (DIM == 3 && SURF) {
1438 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1439 icell = split3d(icell,x);
1440 }
1441 if (DIM < 3 && SURF) {
1442 if (d_cells[icell].nsplit > 1 && d_cells[icell].nsurf >= 0)
1443 icell = split2d(icell,x);
1444 }
1445 } else domain_kk_copy.obj.uncollide_kokkos(outface,x);
1446 } else if (nflag == NPBUNKNOWN) {
1447 icell = -1;
1448 domain_kk_copy.obj.uncollide_kokkos(outface,x);
1449 }
1450
1451 } else if (bflag == SURFACE) {
1452 if (ipart == NULL) {
1453 particle_i.flag = PDISCARD;
1454 break;
1455 } else if (jpart) {
1456 //jpart->flag = PSURF;
1457 //jpart->dtremain = dtremain;
1458 //jpart->weight = particle_i.weight;
1459 //pstop++;
1460 }
1461 Kokkos::atomic_fetch_add(&d_nboundary_one(),1);
1462 Kokkos::atomic_fetch_add(&d_ntouch_one(),-1); // decrement here since will increment below
1463 } else {
1464 if (ATOMIC_REDUCTION == 1) {
1465 Kokkos::atomic_fetch_add(&d_nboundary_one(),1);
1466 Kokkos::atomic_fetch_add(&d_ntouch_one(),-1); // decrement here since will increment below
1467 } else if (ATOMIC_REDUCTION == 0) {
1468 d_nboundary_one()++;
1469 d_ntouch_one()--; // decrement here since will increment below
1470 } else {
1471 reduce.nboundary_one++;
1472 reduce.ntouch_one--; // decrement here since will increment below
1473 }
1474 }
1475 }
1476
1477 // neighbor cell is unknown
1478 // reset icell to original icell which must be a ghost cell
1479 // exit with particle flag = PEXIT, so receiver can identify neighbor
1480
1481 if (icell < 0) {
1482 icell = icell_original;
1483 particle_i.flag = PEXIT;
1484 particle_i.dtremain = dtremain;
1485 d_entryexit() = 1;
1486 break;
1487 }
1488
1489 // if nsurf < 0, new cell is EMPTY ghost
1490 // exit with particle flag = PENTRY, so receiver can continue move
1491
1492 if (d_cells[icell].nsurf < 0) {
1493 particle_i.flag = PENTRY;
1494 particle_i.dtremain = dtremain;
1495 d_entryexit() = 1;
1496 break;
1497 }
1498
1499 // move particle into new grid cell for next stage of move
1500
1501 lo = d_cells[icell].lo;
1502 hi = d_cells[icell].hi;
1503 neigh = d_cells[icell].neigh;
1504 nmask = d_cells[icell].nmask;
1505 if (ATOMIC_REDUCTION == 1)
1506 Kokkos::atomic_fetch_add(&d_ntouch_one(),1);
1507 else if (ATOMIC_REDUCTION == 0)
1508 d_ntouch_one()++;
1509 else
1510 reduce.ntouch_one++;
1511 }
1512
1513 // END of while loop over advection of single particle
1514
1515 #ifdef MOVE_DEBUG
1516 if (ntimestep == MOVE_DEBUG_STEP &&
1517 (MOVE_DEBUG_ID == d_particles[i].id ||
1518 (me == MOVE_DEBUG_PROC && i == MOVE_DEBUG_INDEX)))
1519 printf("MOVE DONE %d %d %d: %g %g %g: DTR %g\n",
1520 MOVE_DEBUG_INDEX,d_particles[i].flag,icell,
1521 x[0],x[1],x[2],dtremain);
1522 #endif
1523
1524 // move is complete, or as much as can be done on this proc
1525 // update particle's grid cell
1526 // if particle flag set, add particle to migrate list
1527 // if discarding, migration will delete particle
1528
1529 particle_i.icell = icell;
1530
1531 if (particle_i.flag != PKEEP) {
1532 int index;
1533 if (ATOMIC_REDUCTION == 0) {
1534 index = d_nmigrate();
1535 d_nmigrate()++;
1536 } else {
1537 index = Kokkos::atomic_fetch_add(&d_nmigrate(),1);
1538 }
1539 k_mlist.d_view[index] = i;
1540 if (particle_i.flag != PDISCARD) {
1541 if (d_cells[icell].proc == me && !d_error_flag()) {
1542 d_error_flag() = 1;
1543 return;
1544 }
1545 if (ATOMIC_REDUCTION == 1)
1546 Kokkos::atomic_fetch_add(&d_ncomm_one(),1);
1547 else if (ATOMIC_REDUCTION == 0)
1548 d_ncomm_one()++;
1549 else
1550 reduce.ncomm_one++;
1551 }
1552 }
1553 } // end of Kokkos parallel_reduce
1554
1555 /* ----------------------------------------------------------------------
1556 particle is entering split parent icell at x
1557 determine which split child cell it is in
1558 return index of sub-cell in ChildCell
1559 ------------------------------------------------------------------------- */
1560
1561 KOKKOS_INLINE_FUNCTION
split3d(int icell,double * x) const1562 int UpdateKokkos::split3d(int icell, double *x) const
1563 {
1564 int m,cflag,isurf,hitflag,side,minsurfindex;
1565 double param,minparam;
1566 double xc[3];
1567 Surf::Tri *tri;
1568
1569 // check for collisions with lines in cell
1570 // find 1st surface hit via minparam
1571 // only consider tris that are mapped via csplits to a split cell
1572 // unmapped tris only touch cell surf at xnew
1573 // another mapped tri should include same xnew
1574 // NOTE: these next 2 lines do not seem correct compared to code
1575 // not considered a collision if particles starts on surf, moving out
1576 // not considered a collision if 2 params are tied and one is INSIDE surf
1577
1578 int nsurf = d_cells[icell].nsurf;
1579 int isplit = d_cells[icell].isplit;
1580 double *xnew = d_sinfo[isplit].xsplit;
1581
1582 cflag = 0;
1583 minparam = 2.0;
1584
1585 auto csplits_begin = d_csplits.row_map(isplit);
1586 auto csurfs_begin = d_csurfs.row_map(icell);
1587 for (m = 0; m < nsurf; m++) {
1588 if (d_csplits.entries(csplits_begin + m) < 0) continue;
1589 isurf = d_csurfs.entries(csurfs_begin + m);
1590 tri = &d_tris[isurf];
1591 hitflag = GeometryKokkos::
1592 line_tri_intersect(x,xnew,
1593 tri->p1,tri->p2,tri->p3,
1594 tri->norm,xc,param,side);
1595
1596 if (hitflag && side != INSIDE && param < minparam) {
1597 cflag = 1;
1598 minparam = param;
1599 minsurfindex = m;
1600 }
1601 }
1602
1603 auto csubs_begin = d_csubs.row_map(isplit);
1604 if (!cflag) return d_csubs.entries(csubs_begin + d_sinfo[isplit].xsub);
1605 int index = d_csplits.entries(csplits_begin + minsurfindex);
1606 return d_csubs.entries(csubs_begin + index);
1607 }
1608
1609 /* ----------------------------------------------------------------------
1610 particle is entering split ICELL at X
1611 determine which split sub-cell it is in
1612 return index of sub-cell in ChildCell
1613 ------------------------------------------------------------------------- */
1614
1615 KOKKOS_INLINE_FUNCTION
split2d(int icell,double * x) const1616 int UpdateKokkos::split2d(int icell, double *x) const
1617 {
1618 int m,cflag,isurf,hitflag,side,minsurfindex;
1619 double param,minparam;
1620 double xc[3];
1621 Surf::Line *line;
1622
1623 // check for collisions with lines in cell
1624 // find 1st surface hit via minparam
1625 // only consider lines that are mapped via csplits to a split cell
1626 // unmapped lines only touch cell surf at xnew
1627 // another mapped line should include same xnew
1628 // NOTE: these next 2 lines do not seem correct compared to code
1629 // not considered a collision if particle starts on surf, moving out
1630 // not considered a collision if 2 params are tied and one is INSIDE surf
1631
1632 int nsurf = d_cells[icell].nsurf;
1633 int isplit = d_cells[icell].isplit;
1634 double *xnew = d_sinfo[isplit].xsplit;
1635
1636 cflag = 0;
1637 minparam = 2.0;
1638 auto csplits_begin = d_csplits.row_map(isplit);
1639 auto csurfs_begin = d_csurfs.row_map(icell);
1640 for (m = 0; m < nsurf; m++) {
1641 if (d_csplits.entries(csplits_begin + m) < 0) continue;
1642 isurf = d_csurfs.entries(csurfs_begin + m);
1643 line = &d_lines[isurf];
1644 hitflag = GeometryKokkos::
1645 line_line_intersect(x,xnew,
1646 line->p1,line->p2,line->norm,
1647 xc,param,side);
1648
1649 if (hitflag && side != INSIDE && param < minparam) {
1650 cflag = 1;
1651 minparam = param;
1652 minsurfindex = m;
1653 }
1654 }
1655
1656 auto csubs_begin = d_csubs.row_map(isplit);
1657 if (!cflag) return d_csubs.entries(csubs_begin + d_sinfo[isplit].xsub);
1658 int index = d_csplits.entries(csplits_begin + minsurfindex);
1659 return d_csubs.entries(csubs_begin + index);
1660 }
1661
1662 /* ----------------------------------------------------------------------
1663 set bounce tally flags for current timestep
1664 nsurf_tally = # of computes needing bounce info on this step
1665 clear accumulators in computes that will be invoked this step
1666 ------------------------------------------------------------------------- */
1667
bounce_set(bigint ntimestep)1668 void UpdateKokkos::bounce_set(bigint ntimestep)
1669 {
1670 Update::bounce_set(ntimestep);
1671
1672 int i;
1673
1674 if (nboundary_tally > KOKKOS_MAX_BLIST)
1675 error->all(FLERR,"Kokkos currently only supports two instances of compute boundary");
1676
1677 if (nboundary_tally) {
1678 for (i = 0; i < nboundary_tally; i++) {
1679 ComputeBoundaryKokkos* compute_boundary_kk = (ComputeBoundaryKokkos*)(blist_active[i]);
1680 compute_boundary_kk->pre_boundary_tally();
1681 blist_active_copy[i].copy(compute_boundary_kk);
1682 }
1683 }
1684
1685 if (nsurf_tally > KOKKOS_MAX_SLIST)
1686 error->all(FLERR,"Kokkos currently only supports two instances of compute surface");
1687
1688 if (nsurf_tally) {
1689 for (i = 0; i < nsurf_tally; i++) {
1690 if (strcmp(slist_active[i]->style,"isurf/grid") == 0)
1691 error->all(FLERR,"Kokkos doesn't yet support compute isurf/grid");
1692 ComputeSurfKokkos* compute_surf_kk = (ComputeSurfKokkos*)(slist_active[i]);
1693 compute_surf_kk->pre_surf_tally();
1694 slist_active_copy[i].copy(compute_surf_kk);
1695 }
1696 }
1697 }
1698