XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 "particles.hpp"
8 #include "space_settings.hpp"
9 #include "distribution.hpp"
10 #include "profile.hpp"
11 
12 extern "C" void set_spall_num_and_ptr(int idx, int n_ptl, int n_vecs, VecParticles* ptl);
13 
14 // Used for Cabana slices (getting one particle property at a time)
15 namespace PtlSlice{
16 enum{Ph=0,Ct,Gid};
17 }
18 
19 // Species class
20 template<class Device>
21 class Species{
22  public:
23 
24  int idx;
25  bool is_electron;
26  bool is_adiabatic;
28  double mass;
29  double charge;
30  double charge_eu;
31  double c_m;
32  double c2_2m;
33 
34  bool is_deltaf;
35 
36  int ncycles;
38 
39  int n_ptl;
40  Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN> particles;
41 
42  // Device particles
43  Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN> particles_d;
46 
47  // phase0 (for ion restoration)
48  Cabana::AoSoA<PhaseDataTypes,HostType,VEC_LEN> phase0;
49  Cabana::AoSoA<PhaseDataTypes,Device,VEC_LEN> phase0_d;
51 
52  // For electron restoration
53  Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN> backup_particles;
54  int n_backup_particles; // Number of particles stored in backup_particles (can't be deduced from its size due to buffer)
55 
57 
58  Eq::Profile<Device> eq_temp; // Equilibrium temperature
59  Eq::Profile<Device> eq_den; // Equilibrium density
60  Eq::Profile<Device> eq_flow; // Equilibrium flow
61 
62  Species(int idx_in, bool is_electron_in, bool is_adiabatic_in, bool is_gyrokinetic_in, double mass_in, double charge_in, double charge_eu_in, bool is_deltaf_in,
63  int ncycles_in);
64 
65  Species(NLReader::NamelistReader& nlr, int idx_in);
66 
67  // Default constructor
68  Species() : owns_particles_d(false) {}
69 
70  // Electron or ion default constructor
71  Species(SpeciesType sp_type, int n_ptl)
72  : idx(sp_type==ELECTRON ? 0 : 1),
73  is_electron(sp_type==ELECTRON),
74  mass(is_electron ? 3.344e-30 : PROTON_MASS),
76  charge_eu(is_electron ? -1.0 : 1.0),
77  is_deltaf(true),
78  is_adiabatic(false),
80  ncycles(is_electron ? 70 : 1),
81  c_m(charge/mass),
82  c2_2m(0.5*charge*charge/mass),
83  n_ptl(n_ptl),
84  backup_particles("backup_particles", 0),
85  particles("particles", add_vec_buffer(n_ptl)),
87  eq_temp(1.0e3,-0.1),
88  eq_den(1.0e19,-0.1),
89  owns_particles_d(false) {}
90 
91  // Special constructor for tests that involve tracking particles in memory
92  // The idea is to use these particles to test functions that reorder particles,
93  // e.g. sort, shift, and cleaning
94  Species(int n_ptl_in)
95  : n_ptl(n_ptl_in),
96  is_adiabatic(false),
97  particles("particles", add_vec_buffer(n_ptl_in))
98  {
99  // Slice particle properties
100  auto ph = Cabana::slice<PtlSlice::Ph>(particles);
101  auto ct = Cabana::slice<PtlSlice::Ct>(particles);
102  auto gid = Cabana::slice<PtlSlice::Gid>(particles);
103 
104  // Assign trackable values
105  for (int i=0;i<n_ptl;i++){
106  // Set GID in order
107  gid(i) = i+1; // 1-indexed
108 
109  // Value of properties is gid + 0.1*(property index)
110  // First particle: (1.0, 1.1, ... 1.8)
111  // Second particle: (2.0, 2.1, ... 2.8)
112  for (int j=0;j<6;j++) ph(i, j) = gid(i) + (j)*0.1;
113  for (int j=0;j<3;j++) ct(i, j) = gid(i) + (j+6)*0.1;
114  }
115 
116  // Buffer particles: same but with gid = -1
117  if(n_ptl>0){
118  for (int i=n_ptl;i<add_vec_buffer(n_ptl);i++){
119  gid(i) = -1;
120  for (int j=0;j<6;j++) ph(i, j) = gid(i) + (j)*0.1;
121  for (int j=0;j<3;j++) ct(i, j) = gid(i) + (j+6)*0.1;
122  }
123  }
124  }
125 
126  void resize_particles(int new_n_ptl){
127  n_ptl = new_n_ptl;
128  particles.resize(add_vec_buffer(n_ptl));
130  }
131 
132  /* If using CPU-only, then "device" particles are a shallow copy of host particles so that
133  * there is no unnecessary duplication. When "device" particles are resized, Cabana will keep the
134  * original allocation if there is a second reference (i.e. host particles).
135  * To resolve this, we free the host particles here so that there is no second reference
136  * */
137  void unassign(){
138 #ifndef USE_GPU
139  particles = Cabana::AoSoA<ParticleDataTypes,HostType,VEC_LEN>();
140 #endif
141  }
142 
143  /* Resizes device particles if on GPU, or just creates a shallow copy if CPU only
144  * */
146  if(!owns_particles_d) exit_XGC("\nSpecies tried to copy particles to device, but doesn't own the device array.");
147 
148 #ifdef USE_GPU
149  // Resize device particles to match host particles
150  particles_d.resize(particles.size());
151 #else
152  // If kernels are on CPU, do shallow copy
154 #endif
155  // Set pointer to correct device memory
156  ptl = (VecParticles*)(particles_d.data());
157  }
158 
159  /* Copies particles to device - deep copy if using GPU, otherwise shallow copy
160  * Also takes the opportunity to set the buffer particles to realistic values
161  * */
163  if(!owns_particles_d) exit_XGC("\nSpecies tried to copy particles to device, but doesn't own the device array.");
164 
165 #ifdef USE_GPU
166  // Resize device particles to match host particles
167  particles_d.resize(particles.size());
168 
169  // Copy to device
170  Cabana::deep_copy(particles_d, particles);
171 #else
172  // If kernels are on CPU, do shallow copy
174 #endif
175  // Set pointer to correct device memory
176  ptl = (VecParticles*)(particles_d.data());
177 
178  // Copy last particle to fill remainder of trailing vector in AoSoA
180  }
181 
182  /* Copies particles from device - deep copy if using GPU, otherwise no copy is necessary
183  * */
185  if(!owns_particles_d) exit_XGC("\nSpecies tried to copy particles from device, but doesn't own the device array.");
186 
187 #ifdef USE_GPU
188  // Copy particles to host
189  Cabana::deep_copy(particles, particles_d);
190 #else
191  // No operation required if CPU-only
192 #endif
193  }
194 
203  if (n_ptl>0){
204  int last_ptl_index = n_ptl - 1;
205  auto ph = Cabana::slice<PtlSlice::Ph>(particles_d);
206  auto ct = Cabana::slice<PtlSlice::Ct>(particles_d);
207  auto gid = Cabana::slice<PtlSlice::Gid>(particles_d);
208 
209  Kokkos::parallel_for("set_buffer_particles_d", Kokkos::RangePolicy<ExSpace>( n_ptl, add_vec_buffer(n_ptl) ), KOKKOS_LAMBDA( const int i ){
210  // Buffer particles: same as last particle, gid = -1
211  for (int j=0;j<6;j++) ph(i, j) = ph(last_ptl_index, j);
212  for (int j=0;j<3;j++) ct(i, j) = ct(last_ptl_index, j);
213  gid(i) = -1;
214  });
215  }
216  }
217 
226  if (n_ptl>0){
227  int last_ptl_index = n_ptl - 1;
228  auto ph = Cabana::slice<PtlSlice::Ph>(phase0_d);
229 
230  Kokkos::parallel_for("set_buffer_phase0", Kokkos::RangePolicy<ExSpace>( n_ptl, add_vec_buffer(n_ptl) ), KOKKOS_LAMBDA( const int i ){
231  // copy final real particle
232  for (int j=0;j<6;j++) ph(i, j) = ph(last_ptl_index, j);
233  });
234  }
235  }
236 };
237 
238 template<class Device>
240  public:
241 
242  TmpSpecies(int n_ptl_in);
244 
245  void resize(int n_ptl_in);
246 
247  // Particles
248  Cabana::AoSoA<ParticleDataTypes,Device,VEC_LEN> particles;
249  VecParticles* ptl; // Pointer to the particles
250 
251  // phase0 (for ions)
252  Cabana::AoSoA<PhaseDataTypes,Device,VEC_LEN> phase0;
253  VecPhase* ph0; // Pointer to phase0
254 };
255 
256 #include "species.tpp"
257 #endif
Cabana::AoSoA< PhaseDataTypes, HostType, VEC_LEN > phase0
Definition: species.hpp:48
Definition: globals.hpp:67
Definition: species.hpp:16
void unassign()
Definition: species.hpp:137
KOKKOS_INLINE_FUNCTION int divide_and_round_up(int a, int b)
Definition: globals.hpp:101
bool owns_particles_d
Whether the species owns the device particle allocation right now.
Definition: species.hpp:45
void set_spall_num_and_ptr(int idx, int n_ptl, int n_vecs, VecParticles *ptl)
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:56
Cabana::AoSoA< ParticleDataTypes, HostType, VEC_LEN > backup_particles
Copy of particles to be restored for RK2.
Definition: species.hpp:53
bool is_electron
Whether this species is the electrons.
Definition: species.hpp:25
double c2_2m
c2/2m
Definition: species.hpp:32
double c_m
c/m
Definition: species.hpp:31
VecPhase * ph0
Definition: species.hpp:253
Definition: species.hpp:16
Eq::Profile< Device > eq_den
Definition: species.hpp:59
Definition: NamelistReader.hpp:163
VecPhase * ph0
Pointer to device phase0.
Definition: species.hpp:50
int add_vec_buffer(int n_ptl)
Definition: particles.hpp:120
int idx
Index in all_species.
Definition: species.hpp:24
TmpSpecies()
Definition: species.hpp:243
Definition: particles.hpp:38
bool is_gyrokinetic
Whether the species is gyrokinetic or drift kinetic.
Definition: species.hpp:27
Definition: distribution.hpp:18
Definition: particles.hpp:55
int n_ptl
Number of particles.
Definition: species.hpp:39
void set_buffer_phase0_d()
Definition: species.hpp:225
void set_buffer_particles_d()
Definition: species.hpp:202
const double UNIT_CHARGE
Charge of an electron (C)
Definition: globals.hpp:93
double charge_eu
Particle charge in eu.
Definition: species.hpp:30
void resize_particles(int new_n_ptl)
Definition: species.hpp:126
VecParticles * ptl
Pointer to the device particles.
Definition: species.hpp:44
double mass
Particle mass.
Definition: species.hpp:28
Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > particles_d
Particles on device.
Definition: species.hpp:43
void resize(int n_ptl_in)
Definition: species.tpp:145
double charge
Particle charge.
Definition: species.hpp:29
void copy_particles_from_device()
Definition: species.hpp:184
int ncycles_between_sorts
Number of subcycles between sorts.
Definition: species.hpp:37
Cabana::AoSoA< PhaseDataTypes, Device, VEC_LEN > phase0_d
Definition: species.hpp:49
Species(SpeciesType sp_type, int n_ptl)
Definition: species.hpp:71
bool is_deltaf
Whether this species is deltaf.
Definition: species.hpp:34
Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > particles
Definition: species.hpp:248
Definition: species.hpp:16
void exit_XGC(std::string msg)
Definition: globals.hpp:26
bool is_adiabatic
Whether this species is adiabatic.
Definition: species.hpp:26
Cabana::AoSoA< PhaseDataTypes, Device, VEC_LEN > phase0
Definition: species.hpp:252
int n_backup_particles
Definition: species.hpp:54
Eq::Profile< Device > eq_flow
Definition: species.hpp:60
void copy_particles_to_device()
Definition: species.hpp:162
Species(int n_ptl_in)
Definition: species.hpp:94
const double PROTON_MASS
Definition: globals.hpp:96
Definition: species.hpp:21
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:234
Eq::Profile< Device > eq_temp
Definition: species.hpp:58
Definition: species.hpp:239
VecParticles * ptl
Definition: species.hpp:249
Species()
Definition: species.hpp:68
int ncycles
Number of subcycles.
Definition: species.hpp:36
Definition: profile.hpp:65
SpeciesType
Definition: globals.hpp:66
void resize_device_particles()
Definition: species.hpp:145
Cabana::AoSoA< ParticleDataTypes, HostType, VEC_LEN > particles
Particles.
Definition: species.hpp:40