XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends 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_ptrs(double* f0_T_ev_cpp, double* f0_fg_T_ev_cpp, double* f0_inv_grid_vol_cpp, double* f0_grid_vol_cpp,
8  double* f0_grid_vol_vonly_cpp, double* f0_den_cpp,
9  double* f0_flow_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 
14 class Plasma{
15 
16  // Variables for managing the device particle memory allocation
20 
21 
22 #ifdef DELTAF_CONV
23  static constexpr bool reduced_deltaf = true;
24 #else
25  static constexpr bool reduced_deltaf = false;
26 #endif
28 
29  public:
30 
31  bool f0_grid;
32 
33  std::vector<Species<DeviceType>> all_species;
34 
35  // Poloidally decomposed f0 values
36  // These will probably become a member of species later, but the fortran code expects
37  // the species to be an array index
39  View<double*,CLayout,HostType> f0_node_cost;
40 
41  private:
42 
43  /* Contains f0 values that are poloidally decomposed but don't need to be transferred between ranks during load rebalancing
44  * */
46  View<double**,CLayout, HostType> f0_T_ev;
47  View<double**,CLayout, HostType> f0_fg_T_ev;
48  View<double**,CLayout, HostType> f0_fg_vth_inv;
49  View<double**,CLayout, HostType> f0_inv_grid_vol;
50  View<double**,CLayout, HostType> f0_grid_vol;
51  View<double**,CLayout, HostType> f0_grid_vol_vonly;
52  View<double**,CLayout, HostType> f0_den;
53  View<double**,CLayout, HostType> f0_flow;
54 
55 
57 
59  : f0_T_ev("f0_T_ev", nsp, pol_decomp.nnodes),
60  f0_fg_T_ev("f0_fg_T_ev", nsp, pol_decomp.nnodes),
61  f0_fg_vth_inv("f0_fg_vth_inv", nsp, pol_decomp.nnodes),
62  f0_inv_grid_vol("f0_inv_grid_vol", nsp, pol_decomp.nnodes),
63  f0_grid_vol("f0_grid_vol", nsp, pol_decomp.nnodes),
64  f0_grid_vol_vonly("f0_grid_vol_vonly", nsp, pol_decomp.nnodes),
65  f0_den("f0_den", nsp, pol_decomp.nnodes),
66  f0_flow("f0_flow", nsp, pol_decomp.nnodes)
67  {
68 #ifndef NO_FORTRAN_MODULES
69  // Fortran arrays are set to point to these Views
71  f0_fg_T_ev.data(),
72  f0_inv_grid_vol.data(),
73  f0_grid_vol.data(),
74  f0_grid_vol_vonly.data(),
75  f0_den.data(),
76  f0_flow.data());
77 #endif
78  }
79  };
80 
81 
82  public:
83 
85 
86  // Not poloidally decomposed
87  View<double**,CLayout, HostType> f0_delta_n;
88  View<double**,CLayout, HostType> f0_delta_u;
89  View<double**,CLayout, HostType> f0_delta_T;
90 
91  View<double**,CLayout, HostType> f0_den_global;
92  View<double**,CLayout, HostType> f0_temp_global;
93 
94  int nspecies;
96 
97  std::vector<std::string> sp_names;
98 
99  // Constructors
101  : particles_d_has_owner(false),
103  nspecies(0), // Initialize with 0 species
104  n_nonadiabatic_species(0), // Initialize with 0 species
105  f0_grid(true),
106  sp_names{"e", "i", "i2", "i3", "i4", "i5", "i6"}
107  {}
108 
110 
114  };
115 
119  };
120 
121  // Loop over all species
122  template<typename F>
123  inline void for_all_species(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
124  for(int isp = 0; isp<all_species.size(); isp++){
125  manage_particle_ownership(isp, device_ptl_opt);
126  func(all_species[isp]);
127  }
128  }
129 
130  // Loop over all non-adiabatic species
131  template<typename F>
132  inline void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
133  for(int isp = 0; isp<all_species.size(); isp++){
134  if(isp>10) { // ISP ERROR
135  printf("ISP ERROR in for_all_nonadiabatic_sepcies: isp=%d",isp);
136  fflush(stdout);
137  }
138  if(!all_species[isp].is_adiabatic){
139  manage_particle_ownership(isp, device_ptl_opt);
140  func(all_species[isp]);
141  }
142  }
143  }
144 
145  // Loop over electrons
146  template<typename F>
147  inline void for_electrons(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
148  for(int isp = 0; isp<all_species.size(); isp++){
149  if(all_species[isp].is_electron){
150  manage_particle_ownership(isp, device_ptl_opt);
151  func(all_species[isp]);
152  }
153  }
154  }
155 
156  // Loop over ions
157  template<typename F>
158  inline void for_all_ions(F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
159  for(int isp = 0; isp<all_species.size(); isp++){
160  if(!all_species[isp].is_electron){
161  manage_particle_ownership(isp, device_ptl_opt);
162  func(all_species[isp]);
163  }
164  }
165  }
166 
167  // Loop over electrons or ions as specified
168  template<typename F>
169  inline void for_all(ParticleType particle_type, F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
170  for(int isp = 0; isp<all_species.size(); isp++){
171  if((particle_type==Electrons && all_species[isp].is_electron) ||
172  (particle_type==Ions && !all_species[isp].is_electron)){
173  manage_particle_ownership(isp, device_ptl_opt);
174  func(all_species[isp]);
175  }
176  }
177  }
178 
179  // Operate on one species
180  template<typename F>
181  inline void for_one_species(int isp, F func, DevicePtlOpt device_ptl_opt = UseDevicePtl){
182  manage_particle_ownership(isp, device_ptl_opt);
183  func(all_species[isp]);
184  }
185 
186  // Loop over all species, return true if the input condition is met for any of them
187  template<typename F>
188  inline bool true_for_some_species(F func) const{
189  for(int isp = 0; isp<all_species.size(); isp++){
190  if(func(all_species[isp])) return true;
191  }
192  return false;
193  }
194 
195  int largest_n_ptl(bool check_backup){
196  int max_n_ptl = 0;
197 
198  // Loop over all species
200  if(species.is_electron && check_backup){
201  // Check on backup particles instead, if that's what's getting used
202  max_n_ptl = std::max(max_n_ptl,species.n_backup_particles);
203  } else {
204  // For now, need to access fortran object for n_ptl
205  max_n_ptl = std::max(max_n_ptl,species.n_ptl);
206  }
207  }, NoDevicePtl);
208 
209  return max_n_ptl;
210  }
211 
212  /* Deallocates the device particles and resets ownership tracking variables
213  * */
215  for_all_nonadiabatic_species([&](Species<DeviceType>& species){
216  if(species.owns_particles_d){
217  species.particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>();
218  species.owns_particles_d = false;
219  }
220  });
221  particles_d_has_owner = false;
222  }
223 
224  static std::vector<MemoryPrediction> estimate_memory_usage(NLReader::NamelistReader& nlr, const Grid<DeviceType> &grid, const DomainDecomposition<DeviceType>& pol_decomp);
225 
226  /* Reallocate and recalculate arrays that are decomposed but recalculated rather than
227  * transferred by the load balance (f0_T_ev etc) */
228  void update_decomposed_f0_calculations(const DomainDecomposition<DeviceType>& pol_decomp,
229  const Grid<DeviceType>& grid,
231  const VelocityGrid& vgrid);
232 
233  private:
234 
235  /* Allocate and calculate global f0 arrays. Should be part of a constructor */
237  // Allocate for all species (even though most simulations don't need it for adiabatic species)
238  f0_delta_n = View<double**,CLayout, HostType>("f0_delta_n", nspecies, grid.nnode);
239  f0_delta_u = View<double**,CLayout, HostType>("f0_delta_u", nspecies, grid.nnode);
240  f0_delta_T = View<double**,CLayout, HostType>("f0_delta_T", nspecies, grid.nnode);
241  // Set all_species unmanaged views
242  for_all_species([&](Species<DeviceType>& species){
243  set_unmanaged_f0_species_view(f0_delta_n, species.idx, species.f0.delta_n_h);
244  set_unmanaged_f0_species_view(f0_delta_u, species.idx, species.f0.delta_u_h);
245  set_unmanaged_f0_species_view(f0_delta_T, species.idx, species.f0.delta_T_h);
246  }, NoDevicePtl);
247 #ifndef NO_FORTRAN_MODULES
248  // Set fortran pointers
249  f0_set_ptrs(grid.nnode, f0_delta_n.data(), f0_delta_u.data(), f0_delta_T.data());
250 #endif
251 
252  // These are constant for the simulation:
253  f0_den_global = View<double**,CLayout, HostType>("f0_den_global", nspecies, grid.nnode);
254  f0_temp_global = View<double**,CLayout, HostType>("f0_temp_global", nspecies, grid.nnode);
255 #ifndef NO_FORTRAN_MODULES
256  // Fortran arrays are set to point to these Views
257  // The calculations are currently done in Fortran too
258  f0_init_global_arrays(f0_den_global.data(), f0_temp_global.data());
259 #endif
260 
261  for_all_species([&](Species<DeviceType>& species){
262  set_unmanaged_f0_species_view(f0_den_global, species.idx, species.f0.den_global_h);
263  set_unmanaged_f0_species_view(f0_temp_global, species.idx, species.f0.temp_global_h);
264 
265  species.calculate_global_f0_arrays(grid, magnetic_field);
266  }, NoDevicePtl);
267  }
268 
269  // Creates an unmanaged view pointing to the subview specified by isp
270  template<typename T_in, typename T_out>
271  void set_unmanaged_f0_species_view(const T_in& view_in, int isp, T_out& view_out){
272  auto view_in_subview = my_subview(view_in, isp);
273  view_out = T_out(view_in_subview.data(), view_in_subview.layout());
274  }
275 
276  public:
277 
278  void resize_f0_f0g(const DomainDecomposition<DeviceType>& pol_decomp, const VelocityGrid& vgrid){
279  // Create new object rather than resizing since we don't need to preserve data
280  f0_f0g = VGridDistribution<HostType>(n_nonadiabatic_species, vgrid, pol_decomp);
281 #ifndef NO_FORTRAN_MODULES
282  int f0_inode1 = pol_decomp.node_offset + 1; // 1-indexed
283  int f0_inode2 = f0_inode1 + pol_decomp.nnodes - 1; // 1-indexed
284  set_f0_f0g_ptr(f0_inode1, f0_inode2, f0_f0g.data());
285 #endif
286 
287  // Set all_species unmanaged views
288  int f0_species_cnt=0;
289  for_all_nonadiabatic_species([&](Species<DeviceType>& species){
290  set_unmanaged_f0_species_view(f0_f0g.f, f0_species_cnt, species.f0.f0g_h);
291  f0_species_cnt++;
292  }, NoDevicePtl);
293  }
294 
295  private:
296 
298  if(!particles_d_has_owner){
299  // Set new owner and initialize with 0 particles
300  particles_d_owner = isp;
301 
302  all_species[particles_d_owner].particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>("particles_d", 0);
303  all_species[particles_d_owner].owns_particles_d = true;
304  particles_d_has_owner = true;
305  }else{
306  // No-op if species is handing off to itself
307  if(particles_d_owner == isp) return;
308 
309  // Shallow copy to new owner
310  all_species[isp].particles_d = all_species[particles_d_owner].particles_d;
311 
312  // Delete original
313  all_species[particles_d_owner].particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>();
314  all_species[particles_d_owner].owns_particles_d = false;
315 
316  // Set new owner
317  particles_d_owner = isp;
318  all_species[particles_d_owner].owns_particles_d = true;
319  }
320  }
321 
322  void manage_particle_ownership(int isp, DevicePtlOpt device_ptl_opt){
323  if((!all_species[isp].is_adiabatic) && device_ptl_opt==UseDevicePtl){
324  if(species_share_particles_d_ownership){
325  transfer_particles_d_ownership(isp);
326  }else{
327  if(!all_species[isp].owns_particles_d){
328  all_species[isp].particles_d = Cabana::AoSoA<ParticleDataTypes,DeviceType,VEC_LEN>("particles_d", 0);
329  all_species[isp].owns_particles_d = true;
330  }
331  }
332 
333  // Resize device particles if they are used
334  all_species[isp].resize_device_particles();
335  }
336  }
337 
338  public:
339  void validate_f0_checkpoint_file_dims(const XGC_IO_Stream& stream, const Grid<DeviceType>& grid, const VelocityGrid& vgrid);
340  void write_checkpoint_files(const Grid<DeviceType>& grid, const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream);
341  void read_checkpoint_files(const Grid<DeviceType>& grid, const VelocityGrid& vgrid, const DomainDecomposition<DeviceType>& pol_decomp, const XGC_IO_Stream& stream, const XGC_IO_Stream& f0_stream, bool n_ranks_is_same, int version);
342 
344  double main_ion_mass = all_species[MAIN_ION].mass;
345  double axis_length = TWOPI * magnetic_field.equil.axis.r;
346  double characteristic_velocity = sqrt(2*main_ion_characteristic_energy/main_ion_mass);
347  return axis_length / characteristic_velocity;
348  }
349 };
350 
351 #endif
void calculate_global_f0_arrays(const Grid< DeviceType > &grid, const MagneticField< DeviceType > &magnetic_field)
Definition: species.cpp:499
void for_all_ions(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:158
bool owns_particles_d
Whether the species owns the device particle allocation right now.
Definition: species.hpp:108
double * data() const
Definition: vgrid_distribution.hpp:73
void init_global_f0_arrays(const Grid< DeviceType > &grid, const MagneticField< DeviceType > &magnetic_field)
Definition: plasma.hpp:236
View< double **, CLayout, HostType > f0_den_global
Equilibrium density at vertices.
Definition: plasma.hpp:91
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:127
void manage_particle_ownership(int isp, DevicePtlOpt device_ptl_opt)
Definition: plasma.hpp:322
bool is_electron
Whether this species is the electrons.
Definition: species.hpp:79
void f0_set_ptrs(int nnode, double *f0_delta_n_cpp, double *f0_delta_u_cpp, double *f0_delta_T_cpp)
static constexpr bool reduced_deltaf
Equivalent to the preprocessor flag for now.
Definition: plasma.hpp:25
ParticleType
Definition: plasma.hpp:116
void for_one_species(int isp, F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:181
View< double **, CLayout, HostType > f0_temp_global
Equilibrium temperature at vertices.
Definition: plasma.hpp:92
Definition: velocity_grid.hpp:8
bool particles_d_has_owner
Whether a species owns the device particles allocation.
Definition: plasma.hpp:18
void deallocate_device_ptl()
Definition: plasma.hpp:214
Definition: plasma.hpp:113
View< double **, CLayout, HostType > f0_fg_vth_inv
thermal velocity from f_grid norm. temp. at nodes (to avoid sqrt ops)
Definition: plasma.hpp:48
Definition: NamelistReader.hpp:193
View< double **, CLayout, HostType > f0_fg_T_ev
f_Grid normalization temperature at nodes
Definition: plasma.hpp:47
Definition: magnetic_field.hpp:13
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:84
int idx
Index in all_species.
Definition: species.hpp:78
VGridDistribution< HostType > f0_f0g
Definition: plasma.hpp:38
bool f0_grid
Definition: plasma.hpp:31
View< double **, CLayout, HostType > f0_flow
Equilibrium flow at nodes.
Definition: plasma.hpp:53
void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:132
View< double **, CLayout, HostType > f0_grid_vol_vonly
Grid volume (v only) at nodes.
Definition: plasma.hpp:51
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:45
View< double **, CLayout, HostType > f0_den
Equilibrium density at nodes.
Definition: plasma.hpp:52
bool default_residence_option()
Definition: species.hpp:35
int n_ptl
Number of particles.
Definition: species.hpp:102
int node_offset
Offset of first mesh node belonging to this MPI rank.
Definition: domain_decomposition.hpp:59
std::vector< Species< DeviceType > > all_species
Every particle species in the simulation.
Definition: plasma.hpp:33
void resize_f0_f0g(const DomainDecomposition< DeviceType > &pol_decomp, const VelocityGrid &vgrid)
Definition: plasma.hpp:278
int nnodes
Number of nodes belonging to this MPI rank.
Definition: domain_decomposition.hpp:60
Definition: plasma.hpp:117
Cabana::AoSoA< ParticleDataTypes, Device, VEC_LEN > particles_d
Particles on device.
Definition: species.hpp:106
View< double **, CLayout, HostType > f0_T_ev
Equilibrium temperature at nodes.
Definition: plasma.hpp:46
bool true_for_some_species(F func) const
Definition: plasma.hpp:188
void set_unmanaged_f0_species_view(const T_in &view_in, int isp, T_out &view_out)
Definition: plasma.hpp:271
int particles_d_owner
Which species, if any, owns the device particles allocation.
Definition: plasma.hpp:19
void f0_init_decomposed_ptrs(double *f0_T_ev_cpp, double *f0_fg_T_ev_cpp, double *f0_inv_grid_vol_cpp, double *f0_grid_vol_cpp, double *f0_grid_vol_vonly_cpp, double *f0_den_cpp, double *f0_flow_cpp)
double get_main_ion_toroidal_transit_time(const MagneticField< DeviceType > &magnetic_field) const
Definition: plasma.hpp:343
void f0_init_global_arrays(double *f0_den_global, double *f0_temp_global)
int largest_n_ptl(bool check_backup)
Definition: plasma.hpp:195
DevicePtlOpt
Definition: plasma.hpp:111
View< double **, CLayout, HostType > f0_delta_T
Flux-surface averaged change of temperature.
Definition: plasma.hpp:89
int nspecies
Number of species including electrons.
Definition: plasma.hpp:94
void for_all_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:123
int n_nonadiabatic_species
Number of nonadiabatic species.
Definition: plasma.hpp:95
View< double **, CLayout, HostType > f0_inv_grid_vol
Inverse grid volume at nodes.
Definition: plasma.hpp:49
void set_f0_f0g_ptr(int f0_inode1, int f0_inode2, double *f0_f0g_loc)
Definition: xgc_io.hpp:24
Definition: plasma.hpp:112
DecomposedRecalculableF0Arrays()
Definition: plasma.hpp:56
double r
Definition: grid_structs.hpp:29
Definition: globals.hpp:85
Kokkos::View< T *, Kokkos::LayoutRight, Device > my_subview(const Kokkos::View< T ****, Kokkos::LayoutRight, Device > &view, int i, int j, int k)
Definition: my_subview.hpp:8
Definition: magnetic_field.F90:1
int n_backup_particles
Definition: species.hpp:124
void for_electrons(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:147
View< double **, CLayout, HostType > f0_delta_n
Flux-surface averaged change of density.
Definition: plasma.hpp:87
Definition: plasma.hpp:14
Definition: plasma.hpp:118
void for_all(ParticleType particle_type, F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:169
bool species_share_particles_d_ownership
Whether to use the device particles sharing scheme.
Definition: plasma.hpp:17
View< double *, CLayout, HostType > f0_node_cost
Definition: plasma.hpp:39
Plasma()
Definition: plasma.hpp:100
Definition: species.hpp:75
void transfer_particles_d_ownership(int isp)
Definition: plasma.hpp:297
View< double **, CLayout, HostType > f0_grid_vol
Grid volume at nodes.
Definition: plasma.hpp:50
View< double **, CLayout, HostType > f0_delta_u
Flux-surface averaged change of parallel flow.
Definition: plasma.hpp:88
std::vector< std::string > sp_names
Definition: plasma.hpp:97
double main_ion_characteristic_energy
Definition: plasma.hpp:27
int nnode
Number of grid nodes.
Definition: grid.hpp:154
RZPair axis
Definition: equil.hpp:86
DecomposedRecalculableF0Arrays(int nsp, const DomainDecomposition< DeviceType > &pol_decomp)
Definition: plasma.hpp:58
constexpr double TWOPI
Definition: constants.hpp:9