XGC1
species.hpp
Go to the documentation of this file.
1 #ifndef SPECIES_HPP
2 #define SPECIES_HPP
3 #include <Cabana_AoSoA.hpp>
4 #include <Cabana_DeepCopy.hpp>
5 #include <Kokkos_Core.hpp>
6 #include "NamelistReader.hpp"
7 #include "timer_macro.hpp"
8 #include "magnetic_field.hpp"
9 #include "grid.hpp"
10 #include "domain_decomposition.hpp"
11 #include "particles.hpp"
12 #include "space_settings.hpp"
13 #include "distribution.hpp"
14 #include "profile.hpp"
15 #include "gyro_avg_mat.hpp"
17 #include "basic_physics.hpp"
18 #include "memory_prediction.hpp"
19 #include "xgc_io.hpp"
20 
21 extern "C" void set_spall_num_and_ptr(int idx, int n_ptl, int n_vecs, VecParticles* ptl);
22 extern "C" void set_min_max_num(int isp, int n_ptl);
23 extern "C" void adjust_n_ptl_for_core_ptl(int* n_ptl);
24 
29  return false;
30 }
31 
36  return false;
37 }
38 
39 // Used for Cabana slices (getting one particle property at a time)
40 namespace PtlSlice{
41 #ifdef ESC_PTL
42 enum{Ph=0,Ct,Gid,Flag};
43 #else
44 enum{Ph=0,Ct,Gid};
45 #endif
46 }
47 
48 struct PtlMvmt{
49  // Options for particle location management when looping over particles
50  enum SendOpt{
51  NoSend=0,
54  };
55 
56  enum ReturnOpt{
60  };
61 
64 
66 };
67 
71 };
72 
73 // Species class
74 template<class Device>
75 class Species{
76  public:
77 
78  int idx;
79  bool is_electron;
80  bool is_adiabatic;
83  double mass;
84  double charge;
85  double charge_eu;
86  double c_m;
87  double c2_2m;
88 
96  bool dynamic_f0;
98 
99  int ncycles;
101 
103  int n_ptl;
104  Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN> particles;
105 
106  // Device particles
107  Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN> particles_d;
108 
110 
113 
114  /*** Could be its own class inside species? ***/
118 
119  // phase0 (for ion restoration)
120  Cabana::AoSoA<PhaseDataTypes,HostType,VEC_LEN> phase0;
121  Cabana::AoSoA<PhaseDataTypes,Device,VEC_LEN> phase0_d;
122 
123  // For electron restoration
124  Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN> backup_particles;
125  Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN> backup_particles_d;
126  int n_backup_particles; // Number of particles stored in backup_particles (can't be deduced from its size due to buffer)
127  bool backup_particles_on_device; // Whether to store back-up particles in device memory
128  /****/
129 
130  View<int*,CLayout,Device> tr_save; // Last known triangle position
131 
133 
135 
136  Eq::Profile<Device> eq_temp; // Equilibrium temperature
137  Eq::Profile<Device> eq_den; // Equilibrium density
138  Eq::Profile<Device> eq_flow; // Equilibrium flow
139  int eq_flow_type; // Type of Equilibirum flow
140 
141  Eq::Profile<Device> eq_fg_temp; // Equilibrium temperature
142  Eq::Profile<Device> eq_fg_flow; // Equilibrium flow - not used for now
143  int eq_fg_flow_type; // Type of Equilibirum flow
144 
145  Eq::Profile<Device> eq_mk_temp; // Equilibrium temperature
146  Eq::Profile<Device> eq_mk_den; // Equilibrium density
147  Eq::Profile<Device> eq_mk_flow; // Equilibrium flow
148  int eq_mk_flow_type; // Type of Equilibirum flow
149 
150 
152 
153  /*** Constructors ***/
154 
155  Species(int idx_in, int nonadiabatic_idx_in, bool is_electron_in, bool is_adiabatic_in, KinType kintype_in, double mass_in, double charge_in, double charge_eu_in,
156  int ncycles_in);
157 
158  Species(NLReader::NamelistReader& nlr, const Grid<DeviceType> &grid, const MagneticField<DeviceType> &magnetic_field, const DomainDecomposition<DeviceType>& pol_decomp, int idx_in, int nonadiabatic_idx_in);
159 
160  // Electron or ion default constructor
161  Species(SpeciesType sp_type, int n_ptl)
162  : idx(sp_type==ELECTRON ? 0 : 1),
163  is_electron(sp_type==ELECTRON),
164  mass(is_electron ? 3.344e-30 : PROTON_MASS),
166  charge_eu(is_electron ? -1.0 : 1.0),
168  nonadiabatic_idx(idx), // Since is_adiabatic is false above
170  ncycles(is_electron ? 70 : 1),
171  c_m(charge/mass),
172  c2_2m(0.5*charge*charge/mass),
174  n_ptl(n_ptl),
175  backup_particles("backup_particles", 0),
176  backup_particles_d("backup_particles_d", 0),
178  particles("particles", add_vec_buffer(n_ptl)),
181  eq_temp(1.0e3,-0.1),
182  eq_den(1.0e19,-0.1),
183  eq_flow_type(2),
184  eq_fg_temp(1.0e3,-0.1),
185  eq_fg_flow_type(2),
186  eq_mk_temp(1.0e3,-0.1),
187  eq_mk_den(1.0e19,-0.1),
188  eq_mk_flow_type(2),
193 
194  // Special constructor for tests that involve tracking particles in memory
195  // The idea is to use these particles to test functions that reorder particles,
196  // e.g. sort, shift, and cleaning
197  Species(int n_ptl_in)
198  : n_ptl(n_ptl_in),
199  is_electron(true),
201  particles("particles", add_vec_buffer(n_ptl_in)),
202  tr_save("tr_save", 0),
209  {
210 
211  // Slice particle properties
212  auto ph = Cabana::slice<PtlSlice::Ph>(particles);
213  auto ct = Cabana::slice<PtlSlice::Ct>(particles);
214  auto gid = Cabana::slice<PtlSlice::Gid>(particles);
215 #ifdef ESC_PTL
216  auto flag = Cabana::slice<PtlSlice::Flag>(particles);
217 #endif
218 
219  // Offset gid if using MPI
220 #ifdef USE_MPI
221  long long int gid_offset = n_ptl*SML_COMM_RANK;
222 #else
223  long long int gid_offset = 0;
224 #endif
225 
226  // Assign trackable values
227  for (int i=0;i<n_ptl;i++){
228  // Set GID in order
229  gid(i) = gid_offset + i+1; // 1-indexed
230 
231  // Value of properties is gid + 0.1*(property index)
232  // First particle: (1.0, 1.1, ... 1.8)
233  // Second particle: (2.0, 2.1, ... 2.8)
234  for (int j=0;j<PTL_NPHASE;j++) ph(i, j) = gid(i) + (j)*0.1;
235  for (int j=0;j<PTL_NCONST;j++) ct(i, j) = gid(i) + (j+PTL_NPHASE)*0.1;
236  }
237 
238  // Buffer particles: same but with gid = -1
239  if(n_ptl>0){
240  for (int i=n_ptl;i<add_vec_buffer(n_ptl);i++){
241  gid(i) = -1;
242  for (int j=0;j<PTL_NPHASE;j++) ph(i, j) = gid(i) + (j)*0.1;
243  for (int j=0;j<PTL_NCONST;j++) ct(i, j) = gid(i) + (j+PTL_NPHASE)*0.1;
244  }
245  }
246  }
247 
248  static std::vector<MemoryPrediction> estimate_memory_usage(NLReader::NamelistReader& nlr, const Grid<DeviceType> &grid, const DomainDecomposition<DeviceType>& pol_decomp, int species_idx);
249 
250  static int get_initial_n_ptl(NLReader::NamelistReader& nlr, const Grid<DeviceType> &grid, const DomainDecomposition<DeviceType>& pol_decomp, int species_idx, bool verbose);
251 
252  void resize_particles(int new_n_ptl){
253  n_ptl = new_n_ptl;
254 
255 #ifndef USE_GPU
256  // If CPU-only, particles_d points to the same location as particles. If particles is resized, then Cabana will not deallocate the first allocation
257  // since it is still used by particles_d. So, reset particles_d before resize, and point it back to particles only afterwards
258  particles_d = Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN>();
259 #endif
260 
261  particles.reserve(minimum_ptl_reservation); // Can only raise reservation (no-op if AoSoA is already larger)
262  particles.resize(add_vec_buffer(n_ptl));
263 
264 #ifndef USE_GPU
265  // Point particles_d back to particles after resize
267 #endif
268 
270  }
271 
273 #ifdef USE_GPU
274  particles.reserve(minimum_ptl_reservation); // Can only raise reservation (no-op if AoSoA is already larger)
275  particles.resize(particles_d.size());
276 #else
277  // Point particles back to particles_d after resize
279 #endif
280 
282  }
283 
284  /* If using CPU-only, then "device" particles are a shallow copy of host particles so that
285  * there is no unnecessary duplication. When "device" particles are resized, Cabana will keep the
286  * original allocation if there is a second reference (i.e. host particles).
287  * To resolve this, we free the host particles here so that there is no second reference
288  * */
290  particles = Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN>();
291  }
292 
293  /* Resizes device particles if on GPU, or just creates a shallow copy if CPU only
294  * */
296  if(!owns_particles_d) exit_XGC("\nSpecies tried to resize device particles, but doesn't own the device array.");
297 
298 #ifdef USE_GPU
299  particles_d.reserve(minimum_ptl_reservation); // Can only raise reservation (no-op if AoSoA is already larger)
300  // Resize device particles to match n particles
302 #else
303  // If kernels are on CPU, do shallow copy
305 #endif
306  }
307 
308  /* Resizes device particles
309  * */
310  void resize_device_particles(int new_n_ptl){
311  if(!owns_particles_d) exit_XGC("\nSpecies tried to resize device particles, but doesn't own the device array.");
312 
313  n_ptl = new_n_ptl;
314 
315  particles_d.reserve(minimum_ptl_reservation); // Can only raise reservation (no-op if AoSoA is already larger)
316 
317  // Resize device particles to match host particles
319  }
320 
321  /* Copies particles to device - deep copy if using GPU, otherwise shallow copy
322  * Also takes the opportunity to set the buffer particles to realistic values
323  * */
325  if(!owns_particles_d) exit_XGC("\nSpecies tried to copy particles to device, but doesn't own the device array.");
326 
327 #ifdef USE_GPU
328  // Copy to device
329  Cabana::deep_copy(particles_d, particles);
330 #else
331  // No operation required if CPU-only
332 #endif
333 
334  // Copy last particle to fill remainder of trailing vector in AoSoA
336  }
337 
338  /* Copies particles from device - deep copy if using GPU, otherwise no copy is necessary
339  * */
341  if(!owns_particles_d) exit_XGC("\nSpecies tried to copy particles from device, but doesn't own the device array.");
342 
343 #ifdef USE_GPU
344  // Copy particles to host
345  Cabana::deep_copy(particles, particles_d);
346 #else
347  // No operation required if CPU-only
348 #endif
349  }
350 
351  /* Copies particles to device if they are resident on the device
352  * */
355  }
356 
357  /* Copies particles from device if they are resident on the device
358  * */
361  }
362 
363  /* Copies particles to device if they are NOT resident on the device
364  * */
367  }
368 
369  /* Copies particles from device if they are NOT resident on the device
370  * */
373  }
374 
383  if (n_ptl>0){
384  int last_ptl_index = n_ptl - 1;
385  auto ph = Cabana::slice<PtlSlice::Ph>(particles_d);
386  auto ct = Cabana::slice<PtlSlice::Ct>(particles_d);
387  auto gid = Cabana::slice<PtlSlice::Gid>(particles_d);
388 #ifdef ESC_PTL
389  auto flag = Cabana::slice<PtlSlice::Flag>(particles_d);
390 #endif
391 
392  Kokkos::parallel_for("set_buffer_particles_d", Kokkos::RangePolicy<ExSpace>( n_ptl, add_vec_buffer(n_ptl) ), KOKKOS_LAMBDA( const int i ){
393  // Buffer particles: same as last particle, gid = -1
394  for (int j=0;j<PTL_NPHASE;j++) ph(i, j) = ph(last_ptl_index, j);
395  for (int j=0;j<PTL_NCONST;j++) ct(i, j) = ct(last_ptl_index, j);
396  gid(i) = -1;
397 #ifdef ESC_PTL
398  flag(i) = flag(last_ptl_index);
399 #endif
400  });
401  }
402  }
403 
412  if (n_ptl>0){
413  int last_ptl_index = n_ptl - 1;
414  auto ph = Cabana::slice<PtlSlice::Ph>(phase0_d);
415 
416  Kokkos::parallel_for("set_buffer_phase0", Kokkos::RangePolicy<ExSpace>( n_ptl, add_vec_buffer(n_ptl) ), KOKKOS_LAMBDA( const int i ){
417  // copy final real particle
418  for (int j=0;j<6;j++) ph(i, j) = ph(last_ptl_index, j);
419  });
420  }
421  }
422 
423  // Options for custom launch bounds since kokkos defaults are suboptimal for electron push kernel
424  enum class LaunchBounds{
425  Default,
426  Custom
427  };
428 
434  template<typename F>
435  inline void for_all_particles(const std::string label, F lambda_func) const {
436  Kokkos::RangePolicy<ExSpace> particle_range_policy( 0, p_range<DeviceType>(n_ptl) );
437  Kokkos::parallel_for(label, Opt::require(particle_range_policy, Async), lambda_func);
438  }
439 
440  inline void back_up_SoA(Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN>& backup_SoA, int offset, int n) const{
441  auto ph_b = Cabana::slice<PtlSlice::Ph>(backup_SoA);
442  auto ct_b = Cabana::slice<PtlSlice::Ct>(backup_SoA);
443  auto gid_b = Cabana::slice<PtlSlice::Gid>(backup_SoA);
444 #ifdef ESC_PTL
445  auto flag_b = Cabana::slice<PtlSlice::Flag>(backup_SoA);
446 #endif
447 
448  auto ph = Cabana::slice<PtlSlice::Ph>(particles_d);
449  auto ct = Cabana::slice<PtlSlice::Ct>(particles_d);
450  auto gid = Cabana::slice<PtlSlice::Gid>(particles_d);
451 #ifdef ESC_PTL
452  auto flag = Cabana::slice<PtlSlice::Flag>(particles_d);
453 #endif
454 
455  Kokkos::parallel_for("backup_first_soa", Kokkos::RangePolicy<ExSpace>( 0, n ), KOKKOS_LAMBDA( const int i ){
456  int i_offset = i + offset;
457  // Make backup copy
458  for (int j=0;j<PTL_NPHASE;j++) ph_b(i, j) = ph(i_offset, j);
459  for (int j=0;j<PTL_NCONST;j++) ct_b(i, j) = ct(i_offset, j);
460  gid_b(i) = gid(i_offset);
461 #ifdef ESC_PTL
462  flag_b(i) = flag(i_offset);
463 #endif
464  // Deactivate
465  gid(i_offset) = -1;
466  });
467  }
468 
469  inline void restore_backup_SoA(Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN>& backup_SoA, int offset, int n) const{
470  auto ph_b = Cabana::slice<PtlSlice::Ph>(backup_SoA);
471  auto ct_b = Cabana::slice<PtlSlice::Ct>(backup_SoA);
472  auto gid_b = Cabana::slice<PtlSlice::Gid>(backup_SoA);
473 #ifdef ESC_PTL
474  auto flag_b = Cabana::slice<PtlSlice::Flag>(backup_SoA);
475 #endif
476 
477  auto ph = Cabana::slice<PtlSlice::Ph>(particles_d);
478  auto ct = Cabana::slice<PtlSlice::Ct>(particles_d);
479  auto gid = Cabana::slice<PtlSlice::Gid>(particles_d);
480 #ifdef ESC_PTL
481  auto flag = Cabana::slice<PtlSlice::Flag>(particles_d);
482 #endif
483 
484  Kokkos::parallel_for("backup_first_soa", Kokkos::RangePolicy<ExSpace>( 0, n ), KOKKOS_LAMBDA( const int i ){
485  int i_offset = i + offset;
486  // Restore from backup copy
487  for (int j=0;j<PTL_NPHASE;j++) ph(i_offset, j) = ph_b(i, j);
488  for (int j=0;j<PTL_NCONST;j++) ct(i_offset, j) = ct_b(i, j);
489  gid(i_offset) = gid_b(i);
490 #ifdef ESC_PTL
491  flag(i_offset) = flag_b(i);
492 #endif
493  });
494  }
495 
501  template<typename F>
502  inline void for_particle_range(int begin_idx, int end_idx, const std::string label, F lambda_func) const {
503  if(end_idx <= begin_idx) return; // Return if range is 0 or less
504 
505  // Still need the subset to line up with the AoSoA vector length
506  int first_soa = begin_idx/VEC_LEN;
507  int n_other_ptl_in_first_soa = begin_idx - first_soa*VEC_LEN;
508  bool first_soa_is_partial = (n_other_ptl_in_first_soa>0);
509  int last_soa = (end_idx-1)/VEC_LEN;
510  int n_other_ptl_in_last_soa = (last_soa+1)*VEC_LEN - end_idx;
511  bool last_soa_is_partial = (n_other_ptl_in_last_soa>0);
512 
513 #ifdef USE_GPU
514  int first_item_in_shifted_range = first_soa*VEC_LEN;
515 #else
516  int first_item_in_shifted_range = first_soa;
517 #endif
518 
519  Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN> ptl_first_soa;
520  Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN> ptl_last_soa;
521  if(first_soa_is_partial){
522  // Make a backup of the first SoA and set the particles_d GIDs to -1
523  ptl_first_soa = Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN>("ptl_first_soa", n_other_ptl_in_first_soa);
524  back_up_SoA(ptl_first_soa, first_soa*VEC_LEN, n_other_ptl_in_first_soa);
525  }
526  if(last_soa_is_partial){
527  // Make a backup of the last SoA and set the particles_d GIDs to -1
528  ptl_last_soa = Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN>("ptl_last_soa", n_other_ptl_in_last_soa);
529  back_up_SoA(ptl_last_soa, end_idx, n_other_ptl_in_last_soa);
530  }
531 
532  // Finally, do the parallel_for
533  Kokkos::RangePolicy<ExSpace> particle_range_policy( first_item_in_shifted_range, p_range<DeviceType>(end_idx) );
534  Kokkos::parallel_for(label, Opt::require(particle_range_policy, Async), lambda_func);
535 
536  // Restore from backup
537  if(first_soa_is_partial){
538  restore_backup_SoA(ptl_first_soa, first_soa*VEC_LEN, n_other_ptl_in_first_soa);
539  }
540  if(last_soa_is_partial){
541  restore_backup_SoA(ptl_last_soa, end_idx, n_other_ptl_in_last_soa);
542  }
543  }
544 
552  template<typename F>
553  inline void for_all_particles(const std::string label, F lambda_func,
554  const PtlMvmt mvmt, LaunchBounds launch_bounds=LaunchBounds::Default) {
555  if(!owns_particles_d) exit_XGC("\nSpecies tried to loop over particles on device, but doesn't own the device array.");
556 
557  bool use_streaming = stream_particles;
558 #ifndef USE_STREAMS
559  use_streaming = false; // Just to be safe, turn streams off here
560 #endif
561 
562  bool send_ptl = ( mvmt.send_opt==PtlMvmt::Send ||
564  bool return_ptl = ( mvmt.return_opt==PtlMvmt::Return ||
566 
567  // Don't need to stream if particles are already present and don't need to be returned
568  if((!send_ptl) && (!return_ptl)) use_streaming = false;
569 
570  if(use_streaming){
571 #ifdef USE_STREAMS
572  Streamed::Option stream_option = Streamed::Normal; // Send to device and back
573  if(!send_ptl) stream_option = Streamed::NoSend;
574  if(!return_ptl) stream_option = Streamed::NoReturn;
575 
576  // Execute streaming parallel_for
577  Streamed::parallel_for(label, n_ptl, lambda_func, stream_option, particles, particles_d);
578 #endif
579  }else{
580  if(send_ptl) TIMER("copy_ptl_to_device",copy_particles_to_device() );
581 
582  if (launch_bounds==LaunchBounds::Custom) {
583 #ifdef USE_EPUSH_LAUNCH_BOUNDS
584 # if !defined(PUSH_MAX_THREADS_PER_BLOCK) || !defined(PUSH_MIN_WARPS_PER_EU)
585 # error "USE_EPUSH_LAUNCH_BOUNDS requires PUSH_MAX_THREADS_PER_BLOCK and PUSH_MIN_WARPS_PER_EU to be defined"
586 # endif
587  Kokkos::RangePolicy<ExSpace, Kokkos::LaunchBounds<PUSH_MAX_THREADS_PER_BLOCK, PUSH_MIN_WARPS_PER_EU>>
588  particle_range_policy( 0, p_range<DeviceType>(n_ptl) );
589  Kokkos::parallel_for(label, Opt::require(particle_range_policy, Async), lambda_func);
590 #else
591  exit_XGC("\nERROR: LaunchBounds::Custom specified, but USE_EPUSH_LAUNCH_BOUNDS is not defined\n");
592 #endif
593  } else {
594  Kokkos::RangePolicy<ExSpace>
595  particle_range_policy( 0, p_range<DeviceType>(n_ptl) );
596  Kokkos::parallel_for(label, Opt::require(particle_range_policy, Async), lambda_func);
597  }
598 
599  if(return_ptl) TIMER("copy_ptl_from_device", copy_particles_from_device() );
600  }
601  }
602 
603  KOKKOS_INLINE_FUNCTION VecParticles* ptl() const{
604  return (VecParticles*)(&particles_d.access(0));
605  }
606 
607  KOKKOS_INLINE_FUNCTION VecPhase* ph0() const{
608  return (VecPhase*)(&phase0_d.access(0));
609  }
610 
612  for_all_particles("copy_to_phase0", KOKKOS_LAMBDA( const int idx ){
614  VecParticles* ptl_loc = species.ptl();
615  VecPhase* ph0_loc = species.ph0();
616  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
617  int p_vec = inds.a + i_simd;
618  ph0_loc[inds.s].r[p_vec] = ptl_loc[inds.s].ph.r[p_vec];
619  ph0_loc[inds.s].z[p_vec] = ptl_loc[inds.s].ph.z[p_vec];
620  ph0_loc[inds.s].phi[p_vec] = ptl_loc[inds.s].ph.phi[p_vec];
621  ph0_loc[inds.s].rho[p_vec] = ptl_loc[inds.s].ph.rho[p_vec];
622  ph0_loc[inds.s].w1[p_vec] = ptl_loc[inds.s].ph.w1[p_vec];
623  ph0_loc[inds.s].w2[p_vec] = ptl_loc[inds.s].ph.w2[p_vec];
624  }
625  });
626  }
627 
628  bool phase0_is_stored() const{
630  }
631 
634  phase0_d.resize(phase0.size());
635  Cabana::deep_copy(phase0_d, phase0);
636  }
637  }
638 
642  backup_particles_d.resize(particles_d.size());
643  Cabana::deep_copy(backup_particles_d, particles_d);
644  }else{
645  // If particles are resident on the device, then use the host particle allocation as the backup
646  // Otherwise, resize the backup_particle array
649  }else{
650  backup_particles.resize(particles_d.size());
651  }
652 
653  // Copy particle data to backup
654  Cabana::deep_copy(backup_particles, particles_d);
655  }
657  } else {
658  // For ions, save phase to phase0
659  phase0_d.resize(particles_d.size());
660  copy_to_phase0(*this); // Separate function to avoid implicit copy into lambda
661  }
663  }
664 
667  // If particles are resident on device, don't need to resize host particles since they are already used as the backup.
668  // If not, resize host particle array to fit backup particles
672  }else{
673  // If particles are not resident on device, resize host particle array to fit backup particles
675  }
676 
677  // Resize device particle array to fit backed up particles
679 
681  Cabana::deep_copy(particles_d, backup_particles_d);
682  }else{
683  // Restore particle data from backup
684  Cabana::deep_copy(particles_d, backup_particles);
685 
687  // If the backup particles are simply pointing to the host particles, then
688  // point them to their own 0-sized allocation when finished using them
689  // so that they don't make a copy when host particles get resized
690  backup_particles = Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN>("backup_particles", 0);
691  }
692  }
693  } else {
695 
696  // Fill buffer with realistic ptl data
698  }
699  particles_are_backed_up = false;
700  }
701 
702  // If phase0 is needed (i.e. for ions), copy phase0 back to host and deallocate device copy
705  phase0.resize(phase0_d.size());
706  Cabana::deep_copy(phase0, phase0_d);
707  phase0_d.resize(0);
708  }
709  }
710 
711  // Deallocate phase0
714  phase0.resize(0);
715  }
716  phase0_d.resize(0);
717  }
718 
719  KOKKOS_INLINE_FUNCTION void restore_phase_from_phase0(const AoSoAIndices<Device>& inds, SimdParticles& part_one) const {
720  VecPhase* ph0_loc = ph0();
721  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
722  int p_vec = inds.a + i_simd;
723  part_one.ph.r[i_simd] = ph0_loc[inds.s].r[p_vec];
724  part_one.ph.z[i_simd] = ph0_loc[inds.s].z[p_vec];
725  part_one.ph.phi[i_simd] = ph0_loc[inds.s].phi[p_vec];
726  part_one.ph.rho[i_simd] = ph0_loc[inds.s].rho[p_vec];
727  part_one.ph.w1[i_simd] = ph0_loc[inds.s].w1[p_vec];
728  part_one.ph.w2[i_simd] = ph0_loc[inds.s].w2[p_vec];
729  }
730  }
731 
732  long long int get_total_n_ptl(){
733 #ifdef USE_MPI
734  long long int tmp_n_ptl = n_ptl;
735  long long int out_n_ptl = 0;
736  MPI_Allreduce(&tmp_n_ptl, &out_n_ptl, 1, MPI_LONG_LONG_INT, MPI_SUM, SML_COMM_WORLD);
737  return out_n_ptl;
738 #else
739  return (long long int)(n_ptl);
740 #endif
741  }
742 
744 #ifdef USE_MPI
745  int tmp_n_ptl = n_ptl;
746  int out_n_ptl = 0;
747  MPI_Allreduce(&tmp_n_ptl, &out_n_ptl, 1, MPI_INT, MPI_MAX, SML_COMM_WORLD);
748  return out_n_ptl;
749 #else
750  return n_ptl;
751 #endif
752  }
753 
754  // Gets the gyro_radius of a species based on equilibrium temperature
755  // inode is the LOCAL (poloidally decomposed) grid node index to get temperature
756  // smu_n is the normalized sqrt(mu)
757  // bfield is the magnetic field at inode
758  KOKKOS_INLINE_FUNCTION double get_fg_gyro_radius(int inode, double smu_n, double bfield) const{
759  // Should replace UNIT_CHARGE*charge_eu with charge(?)
760  return smu_n*sqrt(mass*f0.fg_temp_ev(inode)*EV_2_J) / (UNIT_CHARGE*charge_eu*bfield);
761  }
762 
763  // Gets the equilibrium thermal velocity of a species based on f0 temperature
764  // inode is the GLOBAL node index to get temperature
765  KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity(int inode) const{
766  return thermal_velocity(mass, f0.temp_global(inode));
767  }
768 
769  // Gets the equilibrium thermal velocity of a species based on f0 temperature, on device
770  // inode is the local node index to get temperature
771  KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity_lnode(int inode) const{
772  return thermal_velocity(mass, f0.temp_ev(inode));
773  }
774 
775  // Gets the equilibrium thermal velocity of a species based on f0 temperature, on host
776  // inode is the local node index to get temperature
777  KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity_lnode_h(int inode) const{
778  return thermal_velocity(mass, f0.temp_ev_h(inode));
779  }
780 
781  KOKKOS_INLINE_FUNCTION double get_f0_fg_unit_velocity_lnode_h(int inode) const{
782  return thermal_velocity(mass, f0.fg_temp_ev_h(inode));
783  }
784 
785  // Get species velocity
786  KOKKOS_INLINE_FUNCTION void get_particle_velocity_and_nearest_node(const Grid<DeviceType>& grid, const MagneticField<DeviceType>& magnetic_field, const DomainDecomposition<DeviceType>& pol_decomp, SimdParticles& part, Simd<double>& smu, Simd<double>& vp, Simd<int>& nearest_node, Simd<bool>& not_in_triangle, Simd<bool>& not_in_poloidal_domain) const{
787 
788  // This modulo surely doesnt need to be here (at least, should be elsewhere).
789  // Modulo phi coordinate
790  grid.wedge_modulo_phi(part.ph.phi);
791 
793  grid.get_grid_weights(magnetic_field, part.ph.v(), grid_wts0);
794 
795  // Output argument
796  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
797  not_in_triangle[i_simd] = !grid_wts0.is_valid(i_simd);
798  }
799 
800  Simd<double> bmag;
801  magnetic_field.bmag_interpol(part.ph.v(), bmag);
802 
803  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
804  if(!grid_wts0.is_valid(i_simd)) continue;
805 
806  nearest_node[i_simd]=grid_wts0.node[i_simd] - pol_decomp.node_offset;
807  not_in_poloidal_domain[i_simd] = (nearest_node[i_simd]<0 || nearest_node[i_simd]>=pol_decomp.nnodes);
808 
809  double temp_ev_norm = not_in_poloidal_domain[i_simd] ? f0.fg_temp_ev(0) : f0.fg_temp_ev(nearest_node[i_simd]);
810 
811  // get vp and smu
812  const double& B = bmag[i_simd];
813  vp[i_simd] = normalized_v_para(c_m, mass, B, temp_ev_norm, part.ph.rho[i_simd]);
814  smu[i_simd] = normalized_sqrt_mu(B, temp_ev_norm, part.ct.mu[i_simd]);
815  }
816  }
817 
818  // get flow in m/s according to eq_flow_type
819  KOKKOS_INLINE_FUNCTION double eq_flow_ms(const MagneticField<DeviceType> &magnetic_field, double psi_in,double r,double z, double bphi_over_b) const {
820 
821  double flow = eq_flow.value(magnetic_field,psi_in,r,z);
822  flow=(eq_flow_type>=1)? flow*r : flow;
823  flow=(eq_flow_type>=2)? flow*bphi_over_b : flow ;
824 
825  return flow;
826  }
827 
828  // Return tr_save if available, otherwise return -1
829  KOKKOS_INLINE_FUNCTION void get_tr_save(int i_item, Simd<int>& itr) const{
830  int ibase = i_item*SIMD_SIZE;
831  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
832  int i = ibase + i_simd;
833  itr[i_simd] = i<tr_save.extent(0) ? (tr_save(i)+1) : -1; // tr_save is zero-indexed
834  }
835  }
836 
838  const Grid<DeviceType>& grid,
840  const VelocityGrid& vgrid);
841 
844 
845  void write_ptl_checkpoint_files(const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream, std::string sp_name);
846  void write_f0_checkpoint_files(const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream, std::string sp_name);
847  void read_f0_checkpoint_files(const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream, std::string sp_name);
848  void read_ptl_checkpoint_files(const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream, std::string sp_name, bool n_ranks_is_same, int version);
850 
851  long long int get_max_gid() const;
852  void get_ptl_write_total_and_offsets(const DomainDecomposition<DeviceType>& pol_decomp, long long int& inum_total, long long int& ioff) const;
853 };
854 
855 #endif
KOKKOS_INLINE_FUNCTION double thermal_velocity(double mass, double temp_ev)
Definition: basic_physics.hpp:75
KOKKOS_INLINE_FUNCTION double normalized_v_para(double c_m, double mass, double B, double temp_ev, double rho)
Definition: basic_physics.hpp:97
KOKKOS_INLINE_FUNCTION double normalized_sqrt_mu(double B, double temp_ev, double mu)
Definition: basic_physics.hpp:107
int nnodes
Number of nodes belonging to this MPI rank.
Definition: domain_decomposition.hpp:92
int node_offset
Offset of first mesh node belonging to this MPI rank.
Definition: domain_decomposition.hpp:91
Definition: profile.hpp:171
KOKKOS_INLINE_FUNCTION void wedge_modulo_phi(Simd< double > &phi_mod) const
Definition: grid.tpp:100
KOKKOS_INLINE_FUNCTION void get_grid_weights(const MagneticField< Device > &magnetic_field, const SimdVector &v, const Simd< double > &psi, SimdVector2D &xff, SimdGridWeights< Order::One, PIT > &grid_wts) const
Definition: grid.tpp:32
Definition: gyro_avg_mat.hpp:18
Definition: magnetic_field.hpp:12
Definition: NamelistReader.hpp:193
Definition: species.hpp:75
int n_backup_particles
Definition: species.hpp:126
void read_ptl_checkpoint_files(const DomainDecomposition< DeviceType > &pol_decomp, const XGC_IO_Stream &stream, std::string sp_name, bool n_ranks_is_same, int version)
Definition: species.cpp:671
Eq::Profile< Device > eq_fg_temp
Definition: species.hpp:141
void read_f0_checkpoint_files(const DomainDecomposition< DeviceType > &pol_decomp, const XGC_IO_Stream &stream, std::string sp_name)
Definition: species.cpp:766
Eq::Profile< Device > eq_den
Definition: species.hpp:137
double c2_2m
c2/2m
Definition: species.hpp:87
void copy_phase0_to_device_if_not_resident()
Definition: species.hpp:632
void restore_particles_from_backup()
Definition: species.hpp:665
void resize_device_particles(int new_n_ptl)
Definition: species.hpp:310
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays)
Definition: species.hpp:81
bool backup_particles_on_device
Definition: species.hpp:127
int get_max_n_ptl()
Definition: species.hpp:743
double c_m
c/m
Definition: species.hpp:86
long long int get_max_gid() const
Definition: species.cpp:573
Cabana::AoSoA< ParticleDataTypes, DeviceType, VEC_LEN > backup_particles_d
Copy of particles to be restored for RK2.
Definition: species.hpp:125
View< int *, CLayout, Device > tr_save
Definition: species.hpp:130
KOKKOS_INLINE_FUNCTION VecPhase * ph0() const
Definition: species.hpp:607
KOKKOS_INLINE_FUNCTION double get_f0_fg_unit_velocity_lnode_h(int inode) const
Definition: species.hpp:781
int ncycles
Number of subcycles.
Definition: species.hpp:99
Cabana::AoSoA< PhaseDataTypes, HostType, VEC_LEN > phase0
Definition: species.hpp:120
void resize_host_particles_to_match_device()
Definition: species.hpp:272
int collision_grid_index
Which collision grid to use.
Definition: species.hpp:134
bool particles_are_backed_up
Whether particles are currently backed up.
Definition: species.hpp:117
Species(int idx_in, int nonadiabatic_idx_in, bool is_electron_in, bool is_adiabatic_in, KinType kintype_in, double mass_in, double charge_in, double charge_eu_in, int ncycles_in)
Definition: species.cpp:11
int idx
Index in all_species.
Definition: species.hpp:78
int n_ptl
Number of particles.
Definition: species.hpp:103
void back_up_SoA(Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > &backup_SoA, int offset, int n) const
Definition: species.hpp:440
void save_backup_particles()
Definition: species.hpp:639
void set_buffer_particles_d()
Definition: species.hpp:382
KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity(int inode) const
Definition: species.hpp:765
double charge_eu
Particle charge in eu.
Definition: species.hpp:85
Eq::Profile< Device > eq_flow
Definition: species.hpp:138
Eq::Profile< Device > eq_mk_den
Definition: species.hpp:146
Species(SpeciesType sp_type, int n_ptl)
Definition: species.hpp:161
bool particles_resident_on_device
Whether the particles can reside on device.
Definition: species.hpp:111
MarkerType marker_type
Marker type: reduced delta-f, total-f, full-f, or none (placeholder for adiabatic species)
Definition: species.hpp:89
void copy_to_phase0(Species< Device > &species)
Definition: species.hpp:611
double charge
Particle charge.
Definition: species.hpp:84
int eq_flow_type
Definition: species.hpp:139
Eq::Profile< Device > eq_fg_flow
Definition: species.hpp:142
GyroAverageMatrices< Device > gyro_avg_matrices
Definition: species.hpp:151
bool maxwellian_init
whether initial distribution is maxwellian
Definition: species.hpp:97
KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity_lnode_h(int inode) const
Definition: species.hpp:777
int eq_fg_flow_type
Definition: species.hpp:143
KOKKOS_INLINE_FUNCTION void restore_phase_from_phase0(const AoSoAIndices< Device > &inds, SimdParticles &part_one) const
Definition: species.hpp:719
int ncycles_between_sorts
Number of subcycles between sorts.
Definition: species.hpp:100
int minimum_ptl_reservation
The minimum reservation size for particles.
Definition: species.hpp:102
bool owns_particles_d
Whether the species owns the device particle allocation right now.
Definition: species.hpp:109
double mass
Particle mass.
Definition: species.hpp:83
RKRestorationMethod RK_restoration_method
Currently, electrons must use first method and ions must use second.
Definition: species.hpp:115
KOKKOS_INLINE_FUNCTION double eq_flow_ms(const MagneticField< DeviceType > &magnetic_field, double psi_in, double r, double z, double bphi_over_b) const
Definition: species.hpp:819
KOKKOS_INLINE_FUNCTION void get_tr_save(int i_item, Simd< int > &itr) const
Definition: species.hpp:829
KOKKOS_INLINE_FUNCTION VecParticles * ptl() const
Definition: species.hpp:603
Eq::Profile< Device > eq_mk_temp
Definition: species.hpp:145
LaunchBounds
Definition: species.hpp:424
void restore_backup_SoA(Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > &backup_SoA, int offset, int n) const
Definition: species.hpp:469
Cabana::AoSoA< PhaseDataTypes, Device, VEC_LEN > phase0_d
Definition: species.hpp:121
void copy_particles_to_device_if_not_resident()
Definition: species.hpp:365
void copy_particles_from_device_if_not_resident()
Definition: species.hpp:371
bool is_adiabatic
Whether this species is adiabatic.
Definition: species.hpp:80
Cabana::AoSoA< ParticleDataTypes, HostType, VEC_LEN > particles
Particles.
Definition: species.hpp:104
bool stream_particles
Whether to stream particles between host and device if possible.
Definition: species.hpp:112
void for_particle_range(int begin_idx, int end_idx, const std::string label, F lambda_func) const
Definition: species.hpp:502
KOKKOS_INLINE_FUNCTION void get_particle_velocity_and_nearest_node(const Grid< DeviceType > &grid, const MagneticField< DeviceType > &magnetic_field, const DomainDecomposition< DeviceType > &pol_decomp, SimdParticles &part, Simd< double > &smu, Simd< double > &vp, Simd< int > &nearest_node, Simd< bool > &not_in_triangle, Simd< bool > &not_in_poloidal_domain) const
Definition: species.hpp:786
long long int get_total_n_ptl()
Definition: species.hpp:732
static int get_initial_n_ptl(NLReader::NamelistReader &nlr, const Grid< DeviceType > &grid, const DomainDecomposition< DeviceType > &pol_decomp, int species_idx, bool verbose)
Definition: species.cpp:100
void resize_device_particles()
Definition: species.hpp:295
void copy_particles_to_device()
Definition: species.hpp:324
void for_all_particles(const std::string label, F lambda_func) const
Definition: species.hpp:435
KOKKOS_INLINE_FUNCTION double get_f0_eq_thermal_velocity_lnode(int inode) const
Definition: species.hpp:771
void copy_particles_from_device_if_resident()
Definition: species.hpp:359
void update_decomposed_f0_calculations(const DomainDecomposition< DeviceType > &pol_decomp, const Grid< DeviceType > &grid, const MagneticField< DeviceType > &magnetic_field, const VelocityGrid &vgrid)
Definition: species.cpp:434
KOKKOS_INLINE_FUNCTION double get_fg_gyro_radius(int inode, double smu_n, double bfield) const
Definition: species.hpp:758
void clear_backup_phase()
Definition: species.hpp:712
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:132
Eq::Profile< Device > eq_mk_flow
Definition: species.hpp:147
void move_phase0_from_device_if_not_resident()
Definition: species.hpp:703
static std::vector< MemoryPrediction > estimate_memory_usage(NLReader::NamelistReader &nlr, const Grid< DeviceType > &grid, const DomainDecomposition< DeviceType > &pol_decomp, int species_idx)
Definition: species.cpp:45
void set_buffer_phase0_d()
Definition: species.hpp:411
WeightEvoEq weight_evo_eq
Definition: species.hpp:91
bool phase0_is_stored() const
Definition: species.hpp:628
KinType kintype
Whether the species is gyrokinetic or drift kinetic.
Definition: species.hpp:82
FAnalyticShape f_analytic_shape
f_analytic_shape shape: Maxwellian, SlowingDown or None
Definition: species.hpp:90
Eq::Profile< Device > eq_temp
Definition: species.hpp:136
bool is_electron
Whether this species is the electrons.
Definition: species.hpp:79
Cabana::AoSoA< ParticleDataTypes, HostType, VEC_LEN > backup_particles
Copy of particles to be restored for RK2.
Definition: species.hpp:124
void write_f0_checkpoint_files(const DomainDecomposition< DeviceType > &pol_decomp, const XGC_IO_Stream &stream, std::string sp_name)
Definition: species.cpp:730
void initialize_global_f0_arrays(const Grid< DeviceType > &grid, const MagneticField< DeviceType > &magnetic_field)
Definition: species.cpp:500
Species(int n_ptl_in)
Definition: species.hpp:197
void for_all_particles(const std::string label, F lambda_func, const PtlMvmt mvmt, LaunchBounds launch_bounds=LaunchBounds::Default)
Definition: species.hpp:553
void unassign_host_particles()
Definition: species.hpp:289
void copy_particles_from_device()
Definition: species.hpp:340
void copy_particles_to_device_if_resident()
Definition: species.hpp:353
void resize_particles(int new_n_ptl)
Definition: species.hpp:252
bool dynamic_f0
Whether f0 can evolve in time.
Definition: species.hpp:96
int eq_mk_flow_type
Definition: species.hpp:148
Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > particles_d
Particles on device.
Definition: species.hpp:107
void get_ptl_write_total_and_offsets(const DomainDecomposition< DeviceType > &pol_decomp, long long int &inum_total, long long int &ioff) const
Definition: species.cpp:596
void write_ptl_checkpoint_files(const DomainDecomposition< DeviceType > &pol_decomp, const XGC_IO_Stream &stream, std::string sp_name)
Definition: species.cpp:621
void read_initial_distribution(NLReader::NamelistReader &nlr, const DomainDecomposition< DeviceType > &pol_decomp)
Definition: species.cpp:803
Definition: xgc_io.hpp:24
constexpr double EV_2_J
Conversion rate ev to J.
Definition: constants.hpp:5
constexpr double UNIT_CHARGE
Charge of an electron (C)
Definition: constants.hpp:4
constexpr double PROTON_MASS
Definition: constants.hpp:7
void exit_XGC(std::string msg)
Definition: globals.hpp:37
@ PTL_NPHASE
Definition: globals.hpp:212
FAnalyticShape
Definition: globals.hpp:116
SpeciesType
Definition: globals.hpp:83
@ ELECTRON
Definition: globals.hpp:84
KOKKOS_INLINE_FUNCTION int divide_and_round_up(int a, int b)
Definition: globals.hpp:226
WeightEvoEq
Definition: globals.hpp:121
KinType
Definition: globals.hpp:88
@ GyroKin
Definition: globals.hpp:90
@ DriftKin
Definition: globals.hpp:89
@ PTL_NCONST
Definition: globals.hpp:222
MarkerType
Definition: globals.hpp:110
int SML_COMM_RANK
Definition: my_mpi.cpp:5
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
Definition: species.hpp:40
@ Ph
Definition: species.hpp:44
@ Gid
Definition: species.hpp:44
@ Ct
Definition: species.hpp:44
Option
Definition: streamed_parallel_for.hpp:13
@ NoReturn
Definition: streamed_parallel_for.hpp:16
@ NoSend
Definition: streamed_parallel_for.hpp:14
@ Normal
Definition: streamed_parallel_for.hpp:15
void parallel_for(const std::string name, int n_ptl, Function func, Option option, HostAoSoA aosoa_h, DeviceAoSoA aosoa_d)
Definition: streamed_parallel_for.hpp:252
Definition: magnetic_field.F90:1
logical false
Definition: module.F90:102
logical true
Definition: module.F90:102
int p_range< DeviceType >(int num_particle)
Definition: particles.hpp:193
int add_vec_buffer(int n_ptl)
Definition: particles.hpp:200
constexpr static const Kokkos::Experimental::WorkItemProperty::HintLightWeight_t Async
Definition: space_settings.hpp:83
void set_min_max_num(int isp, int n_ptl)
bool default_residence_option()
Definition: species.hpp:35
RKRestorationMethod
Definition: species.hpp:68
@ RestoreOnOriginalRank
Definition: species.hpp:69
@ BringOldPhaseToNewRank
Definition: species.hpp:70
bool default_streaming_option()
Definition: species.hpp:28
void set_spall_num_and_ptr(int idx, int n_ptl, int n_vecs, VecParticles *ptl)
void adjust_n_ptl_for_core_ptl(int *n_ptl)
Definition: particles.hpp:154
int s
The index in the outer array of the AoSoA.
Definition: particles.hpp:155
int a
The index in the inner array of the AoSoA.
Definition: particles.hpp:156
Definition: distribution.hpp:11
Definition: species.hpp:48
SendOpt send_opt
Definition: species.hpp:62
ReturnOpt return_opt
Definition: species.hpp:63
SendOpt
Definition: species.hpp:50
@ NoSend
Definition: species.hpp:51
@ SendIfNotResident
Definition: species.hpp:53
@ Send
Definition: species.hpp:52
PtlMvmt(SendOpt send_opt, ReturnOpt return_opt)
Definition: species.hpp:65
ReturnOpt
Definition: species.hpp:56
@ ReturnIfNotResident
Definition: species.hpp:59
@ Return
Definition: species.hpp:58
@ NoReturn
Definition: species.hpp:57
Simd< double > mu
m*v_perp^2/(2B)
Definition: particles.hpp:52
Definition: grid_weights.hpp:47
Definition: particles.hpp:61
SimdPhase ph
Definition: particles.hpp:62
SimdConstants ct
Definition: particles.hpp:63
Simd< double > w2
(1 - background distribution)/f0
Definition: particles.hpp:23
Simd< double > w1
delta-f weight
Definition: particles.hpp:22
Simd< double > rho
m*v_para/(q*B) - A_para^h/B (should it be plus or minus?)
Definition: particles.hpp:21
Simd< double > z
Cylindrical coordinate Z.
Definition: particles.hpp:19
Simd< double > r
Cylindrical coordinate R (major radial direction)
Definition: particles.hpp:18
Simd< double > phi
Cylindrical coordinate phi (toroidal direction)
Definition: particles.hpp:20
KOKKOS_INLINE_FUNCTION SimdVector & v()
Definition: particles.hpp:39
Definition: particles.hpp:115
VecPhase ph
Definition: particles.hpp:116
Definition: particles.hpp:95
double r[VEC_LEN]
Definition: particles.hpp:96
double w2[VEC_LEN]
Definition: particles.hpp:101
double phi[VEC_LEN]
Definition: particles.hpp:98
double rho[VEC_LEN]
Definition: particles.hpp:99
double w1[VEC_LEN]
Definition: particles.hpp:100
double z[VEC_LEN]
Definition: particles.hpp:97
Definition: velocity_grid.hpp:8
#define TIMER(N, F)
Definition: timer_macro.hpp:24