XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
plasma.hpp
Go to the documentation of this file.
1 #ifndef PLASMA_HPP
2 #define PLASMA_HPP
3 #include "my_subview.hpp"
4 #include "species.hpp"
5 #include "vgrid_distribution.hpp"
6 
7 extern "C" void f0_init_decomposed_arrays(double* f0_T_ev_cpp, double* f0_inv_grid_vol_cpp, double* f0_grid_vol_cpp,
8  double* f0_grid_vol_vonly_cpp, double* f0_n_Ta_cpp, double* f0_den_cpp,
9  double* f0_flow_cpp, double* f0_B_B0_cpp);
10 extern "C" void f0_init_global_arrays(double* f0_den_global, double* f0_temp_global);
11 extern "C" void f0_set_ptrs(int nnode, double* f0_delta_n_cpp, double* f0_delta_u_cpp, double* f0_delta_T_cpp);
12 extern "C" void set_f0_f0g_ptr(int f0_inode1, int f0_inode2, double* f0_f0g_loc);
13 extern "C" void f0_initialize_wrap();
14 
15 class Plasma{
16 
17  // Variables for managing the device particle memory allocation
21 
22  public:
23 
24  std::vector<Species<DeviceType>> all_species;
25 
26  // Poloidally decomposed f0 values
27  // These will probably become a member of species later, but the fortran code expects
28  // the species to be an array index
30 
31  private:
32  /* Contains f0 values that are poloidally decomposed but don't need to be transferred between ranks during load rebalancing
33  * */
35  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_T_ev;
36  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_inv_grid_vol;
37  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_grid_vol;
38  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_grid_vol_vonly;
39  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_n_Ta;
40  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_den;
41  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_flow;
42 
43  // f0_B_B0 is not plasma-related, should be moved elsewhere (e.g. Grid) - ALS
44  Kokkos::View<double*,Kokkos::LayoutRight, HostType> f0_B_B0;
45 
47 
49  : f0_T_ev("f0_T_ev", nsp, pol_decomp.nnodes),
50  f0_inv_grid_vol("f0_inv_grid_vol", nsp, pol_decomp.nnodes),
51  f0_grid_vol("f0_grid_vol", nsp, pol_decomp.nnodes),
52  f0_grid_vol_vonly("f0_grid_vol_vonly", nsp, pol_decomp.nnodes),
53  f0_n_Ta("f0_n_Ta", nsp, pol_decomp.nnodes),
54  f0_den("f0_den", nsp, pol_decomp.nnodes),
55  f0_flow("f0_flow", nsp, pol_decomp.nnodes),
56  f0_B_B0("f0_B_B0", pol_decomp.nnodes)
57  {
58  // Fortran arrays are set to point to these Views
59  // The calculations are currently done in Fortran too
61  f0_inv_grid_vol.data(),
62  f0_grid_vol.data(),
63  f0_grid_vol_vonly.data(),
64  f0_n_Ta.data(),
65  f0_den.data(),
66  f0_flow.data(),
67  f0_B_B0.data());
68  }
69  };
70 
71  public:
72 
74 
75  // Not poloidally decomposed
76  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_delta_n;
77  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_delta_u;
78  Kokkos::View<double**,Kokkos::LayoutRight, HostType> f0_delta_T;
79 
80  View<double**,CLayout, HostType> f0_den_global;
81  View<double**,CLayout, HostType> f0_temp_global;
82 
83  int nspecies;
85 
86  // Constructors
88  : particles_d_has_owner(false),
90  nspecies(0), // Initialize with 0 species
91  n_nonadiabatic_species(0) // Initialize with 0 species
92  {}
93 
94  Plasma(NLReader::NamelistReader& nlr, bool use_f0_grid, bool allocate_for_adiabatic_electrons, const DomainDecomposition<DeviceType>& pol_decomp, const VelocityGrid& vgrid)
95  : particles_d_has_owner(false)
96  {
97  if(nlr.namelist_present("performance_param")){
98  nlr.use_namelist("performance_param");
99  species_share_particles_d_ownership = !(nlr.get<bool>("particles_resident_on_device", default_residence_option()));
100  }else{
102  }
103 #ifndef USE_GPU
105 #endif
106 
107  nlr.use_namelist("ptl_param");
108  nspecies = nlr.get<int>("ptl_nsp",1);
109  nspecies += 1; // Add electrons
110 
111  // Add each species
112  int i_nonadiabatic_species = 0;
113  for(int isp = 0; isp<nspecies; isp++){
114  all_species.push_back(Species<DeviceType>(nlr, isp, i_nonadiabatic_species));
115  if(!all_species[isp].is_adiabatic){
116  i_nonadiabatic_species++;
117  }
118  }
119 
120  // Save count of non-adiabatic species
121  n_nonadiabatic_species = i_nonadiabatic_species;
122 
123  /*** f0 arrays ***/
124  int grid_nnode = pol_decomp.gvid0_pid_h(pol_decomp.gvid0_pid_h.size()-1) - 1;
125  init_global_f0_arrays(grid_nnode, allocate_for_adiabatic_electrons);
126 
127  if(use_f0_grid){
128  resize_f0_f0g(pol_decomp, vgrid);
130 
131  // Reallocate and recalculate arrays that are decomposed but recalculated rather than
132  // transferred by the load balance (f0_T_ev etc)
134  }
135  }
136 
140  };
141 
145  };
146 
147  // Loop over all species
148  template<typename F>
149  inline void for_all_species(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
150  for(int isp = 0; isp<all_species.size(); isp++){
151  manage_particle_ownership(isp, device_ptl_opt);
152  func(all_species[isp]);
153  }
154  }
155 
156  // Loop over all non-adiabatic species
157  template<typename F>
158  inline void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
159  for(int isp = 0; isp<all_species.size(); isp++){
160  if(isp>10) { // ISP ERROR
161  printf("ISP ERROR in for_all_nonadiabatic_sepcies: isp=%d",isp);
162  fflush(stdout);
163  }
164  if(!all_species[isp].is_adiabatic){
165  manage_particle_ownership(isp, device_ptl_opt);
166  func(all_species[isp]);
167  }
168  }
169  }
170 
171  // Loop over electrons
172  template<typename F>
173  inline void for_electrons(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
174  for(int isp = 0; isp<all_species.size(); isp++){
175  if(all_species[isp].is_electron){
176  manage_particle_ownership(isp, device_ptl_opt);
177  func(all_species[isp]);
178  }
179  }
180  }
181 
182  // Loop over ions
183  template<typename F>
184  inline void for_all_ions(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
185  for(int isp = 0; isp<all_species.size(); isp++){
186  if(!all_species[isp].is_electron){
187  manage_particle_ownership(isp, device_ptl_opt);
188  func(all_species[isp]);
189  }
190  }
191  }
192 
193  // Loop over electrons or ions as specified
194  template<typename F>
195  inline void for_all(ParticleType particle_type, F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
196  for(int isp = 0; isp<all_species.size(); isp++){
197  if((particle_type==Electrons && all_species[isp].is_electron) ||
198  (particle_type==Ions && !all_species[isp].is_electron)){
199  manage_particle_ownership(isp, device_ptl_opt);
200  func(all_species[isp]);
201  }
202  }
203  }
204 
205  // Operate on one species
206  template<typename F>
207  inline void for_one_species(int isp, F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
208  manage_particle_ownership(isp, device_ptl_opt);
209  func(all_species[isp]);
210  }
211 
212  int largest_n_ptl(bool check_backup){
213  int max_n_ptl = 0;
214 
215  // Loop over all species
217  if(species.is_electron && check_backup){
218  // Check on backup particles instead, if that's what's getting used
219  max_n_ptl = std::max(max_n_ptl,species.n_backup_particles);
220  } else {
221  // For now, need to access fortran object for n_ptl
222  max_n_ptl = std::max(max_n_ptl,species.n_ptl);
223  }
224  }, NoDevicePtl);
225 
226  return max_n_ptl;
227  }
228 
229  /* Deallocates the device particles and resets ownership tracking variables
230  * */
232  for_all_nonadiabatic_species([&](Species<DeviceType>& species){
233  if(species.owns_particles_d){
234  species.particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>();
235  species.owns_particles_d = false;
236  }
237  });
238  particles_d_has_owner = false;
239  }
240 
241  /* Reallocate and recalculate arrays that are decomposed but recalculated rather than
242  * transferred by the load balance (f0_T_ev etc) */
244  decomposed_recalculable_f0_arrays = DecomposedRecalculableF0Arrays(n_nonadiabatic_species, pol_decomp);
245 
246  // Set all_species unmanaged views
247  int f0_species_cnt=0;
248  for_all_nonadiabatic_species([&](Species<DeviceType>& species){
249  set_unmanaged_f0_species_view(decomposed_recalculable_f0_arrays.f0_T_ev, f0_species_cnt, species.f0.temp_ev_h);
250  f0_species_cnt++;
251  }, NoDevicePtl);
252  }
253 
254  private:
255 
256  /* Allocate and calculate global f0 arrays. Should be part of a constructor */
257  void init_global_f0_arrays(int grid_nnode, bool is_XGCa){
258  // f0_delta_{n,u,T} are sized differently in XGC1 vs XGCa (the latter uses index 0 even if
259  // electrons are adiabatic)
260  int nsp = (is_XGCa ? nspecies : n_nonadiabatic_species);
261 
262  f0_delta_n = View<double**,CLayout, HostType>("f0_delta_n", nsp, grid_nnode);
263  f0_delta_u = View<double**,CLayout, HostType>("f0_delta_u", nsp, grid_nnode);
264  f0_delta_T = View<double**,CLayout, HostType>("f0_delta_T", nsp, grid_nnode);
265  f0_set_ptrs(grid_nnode, f0_delta_n.data(), f0_delta_u.data(), f0_delta_T.data());
266 
267  // These are constant for the simulation:
268  f0_den_global = View<double**,CLayout, HostType>("f0_den_global", n_nonadiabatic_species, grid_nnode);
269  f0_temp_global = View<double**,CLayout, HostType>("f0_temp_global", n_nonadiabatic_species, grid_nnode);
270 
271  // Fortran arrays are set to point to these Views
272  // The calculations are currently done in Fortran too
273  f0_init_global_arrays(f0_den_global.data(), f0_temp_global.data());
274 
275  for_all_nonadiabatic_species([&](Species<DeviceType>& species){
276  set_unmanaged_f0_species_view(f0_den_global, species.nonadiabatic_idx, species.f0.den_global_h);
277  set_unmanaged_f0_species_view(f0_temp_global, species.nonadiabatic_idx, species.f0.temp_global_h);
278  }, NoDevicePtl);
279  }
280 
281  // Creates an unmanaged view pointing to the subview specified by isp
282  template<typename T_in, typename T_out>
283  void set_unmanaged_f0_species_view(const T_in& view_in, int isp, T_out& view_out){
284  auto view_in_subview = my_subview(view_in, isp);
285  view_out = T_out(view_in_subview.data(), view_in_subview.layout());
286  }
287 
288  public:
289 
290  void resize_f0_f0g(const DomainDecomposition<DeviceType>& pol_decomp, const VelocityGrid& vgrid){
291  // Create new object rather than resizing since we don't need to preserve data
292  f0_f0g = VGridDistribution<HostType>(n_nonadiabatic_species, vgrid, pol_decomp);
293 
294  int f0_inode1 = pol_decomp.first_node; // 1-indexed
295  int f0_inode2 = f0_inode1 + pol_decomp.nnodes - 1; // 1-indexed
296  set_f0_f0g_ptr(f0_inode1, f0_inode2, f0_f0g.data());
297 
298  // Set all_species unmanaged views
299  int f0_species_cnt=0;
300  for_all_nonadiabatic_species([&](Species<DeviceType>& species){
301  set_unmanaged_f0_species_view(f0_f0g.f, f0_species_cnt, species.f0.f0g_h);
302  f0_species_cnt++;
303  }, NoDevicePtl);
304  }
305 
306  private:
307 
309  if(!particles_d_has_owner){
310  // Set new owner and initialize with 0 particles
311  particles_d_owner = isp;
312 
313  all_species[particles_d_owner].particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>("particles_d", 0);
314  all_species[particles_d_owner].owns_particles_d = true;
315  particles_d_has_owner = true;
316  }else{
317  // No-op if species is handing off to itself
318  if(particles_d_owner == isp) return;
319 
320  // Shallow copy to new owner
321  all_species[isp].particles_d = all_species[particles_d_owner].particles_d;
322 
323  // Delete original
324  all_species[particles_d_owner].particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>();
325  all_species[particles_d_owner].owns_particles_d = false;
326 
327  // Set new owner
328  particles_d_owner = isp;
329  all_species[particles_d_owner].owns_particles_d = true;
330  }
331  }
332 
333  void manage_particle_ownership(int isp, DevicePtlOpt device_ptl_opt){
334  if((!all_species[isp].is_adiabatic) && device_ptl_opt==UseDevicePtl){
335  if(species_share_particles_d_ownership){
336  transfer_particles_d_ownership(isp);
337  }else{
338  if(!all_species[isp].owns_particles_d){
339  all_species[isp].particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>("particles_d", 0);
340  all_species[isp].owns_particles_d = true;
341  }
342  }
343 
344  // Resize device particles if they are used
345  all_species[isp].resize_device_particles();
346  }
347  }
348 
349 };
350 
351 #endif
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_n_Ta
Equilibrium n_Ta at nodes.
Definition: plasma.hpp:39
void for_all_ions(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:184
Kokkos::View< double *, Kokkos::LayoutRight, HostType > f0_B_B0
Bfield normalized to B0 at nodes.
Definition: plasma.hpp:44
bool owns_particles_d
Whether the species owns the device particle allocation right now.
Definition: species.hpp:92
View< double **, CLayout, HostType > f0_den_global
Equilibrium density at vertices.
Definition: plasma.hpp:80
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:111
void manage_particle_ownership(int isp, DevicePtlOpt device_ptl_opt)
Definition: plasma.hpp:333
bool is_electron
Whether this species is the electrons.
Definition: species.hpp:70
void f0_set_ptrs(int nnode, double *f0_delta_n_cpp, double *f0_delta_u_cpp, double *f0_delta_T_cpp)
ParticleType
Definition: plasma.hpp:142
void for_one_species(int isp, F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:207
void update_decomposed_f0_calculations(const DomainDecomposition< DeviceType > &pol_decomp)
Definition: plasma.hpp:243
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:353
View< double **, CLayout, HostType > f0_temp_global
Equilibrium temperature at vertices.
Definition: plasma.hpp:81
Definition: velocity_grid.hpp:7
bool particles_d_has_owner
Whether a species owns the device particles allocation.
Definition: plasma.hpp:19
void deallocate_device_ptl()
Definition: plasma.hpp:231
Definition: plasma.hpp:139
Definition: NamelistReader.hpp:163
DecomposedRecalculableF0Arrays decomposed_recalculable_f0_arrays
Contains f0 values that are poloidally decomposed but don&#39;t need to be transferred between ranks duri...
Definition: plasma.hpp:73
VGridDistribution< HostType > f0_f0g
Definition: plasma.hpp:29
void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:158
void f0_initialize_wrap()
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays) ...
Definition: species.hpp:72
bool default_residence_option()
Definition: species.hpp:26
int n_ptl
Number of particles.
Definition: species.hpp:86
Kokkos::View< T *, Kokkos::LayoutRight, Device > my_subview(const Kokkos::View< T ***, Kokkos::LayoutRight, Device > &view, int i, int j)
Definition: my_subview.hpp:8
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_delta_n
Flux-surface averaged change of density.
Definition: plasma.hpp:76
std::vector< Species< DeviceType > > all_species
Every particle species in the simulation.
Definition: plasma.hpp:24
void resize_f0_f0g(const DomainDecomposition< DeviceType > &pol_decomp, const VelocityGrid &vgrid)
Definition: plasma.hpp:290
int nnodes
Number of nodes belonging to this MPI rank.
Definition: domain_decomposition.hpp:39
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_flow
Equilibrium flow at nodes.
Definition: plasma.hpp:41
Definition: plasma.hpp:143
Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > particles_d
Particles on device.
Definition: species.hpp:90
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_inv_grid_vol
Inverse grid volume at nodes.
Definition: plasma.hpp:36
void set_unmanaged_f0_species_view(const T_in &view_in, int isp, T_out &view_out)
Definition: plasma.hpp:283
int particles_d_owner
Which species, if any, owns the device particles allocation.
Definition: plasma.hpp:20
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:322
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_grid_vol_vonly
Grid volume (v only) at nodes.
Definition: plasma.hpp:38
void f0_init_global_arrays(double *f0_den_global, double *f0_temp_global)
int largest_n_ptl(bool check_backup)
Definition: plasma.hpp:212
DevicePtlOpt
Definition: plasma.hpp:137
int nspecies
Number of species including electrons.
Definition: plasma.hpp:83
void for_all_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:149
int n_nonadiabatic_species
Number of nonadiabatic species.
Definition: plasma.hpp:84
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_T_ev
Equilibrium temperature at nodes.
Definition: plasma.hpp:35
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_grid_vol
Grid volume at nodes.
Definition: plasma.hpp:37
void set_f0_f0g_ptr(int f0_inode1, int f0_inode2, double *f0_f0g_loc)
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_delta_u
Flux-surface averaged change of parallel flow.
Definition: plasma.hpp:77
bool namelist_present(const string &namelist)
Definition: NamelistReader.hpp:312
void init_global_f0_arrays(int grid_nnode, bool is_XGCa)
Definition: plasma.hpp:257
Definition: plasma.hpp:138
DecomposedRecalculableF0Arrays()
Definition: plasma.hpp:46
void f0_init_decomposed_arrays(double *f0_T_ev_cpp, double *f0_inv_grid_vol_cpp, double *f0_grid_vol_cpp, double *f0_grid_vol_vonly_cpp, double *f0_n_Ta_cpp, double *f0_den_cpp, double *f0_flow_cpp, double *f0_B_B0_cpp)
int n_backup_particles
Definition: species.hpp:108
void for_electrons(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:173
double * data()
Definition: vgrid_distribution.hpp:39
Definition: plasma.hpp:15
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_den
Equilibrium density at nodes.
Definition: plasma.hpp:40
Definition: plasma.hpp:144
void for_all(ParticleType particle_type, F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:195
bool species_share_particles_d_ownership
Whether to use the device particles sharing scheme.
Definition: plasma.hpp:18
Plasma()
Definition: plasma.hpp:87
Definition: species.hpp:66
void transfer_particles_d_ownership(int isp)
Definition: plasma.hpp:308
Kokkos::View< double **, Kokkos::LayoutRight, HostType > f0_delta_T
Flux-surface averaged change of temperature.
Definition: plasma.hpp:78
int first_node
First mesh node belonging to this MPI rank.
Definition: domain_decomposition.hpp:38
DecomposedRecalculableF0Arrays(int nsp, const DomainDecomposition< DeviceType > &pol_decomp)
Definition: plasma.hpp:48
Plasma(NLReader::NamelistReader &nlr, bool use_f0_grid, bool allocate_for_adiabatic_electrons, const DomainDecomposition< DeviceType > &pol_decomp, const VelocityGrid &vgrid)
Definition: plasma.hpp:94
Kokkos::View< int *, Kokkos::LayoutRight, HostType > gvid0_pid_h
Which processors get which vertices (host)
Definition: domain_decomposition.hpp:42