XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
col_species.hpp
Go to the documentation of this file.
1 #ifndef COL_SPECIES_HPP
2 #define COL_SPECIES_HPP
3 
4 #include <vector>
5 #include <Kokkos_Core.hpp>
6 #include <Kokkos_DualView.hpp>
7 
8 #include "vgrid_distribution.hpp"
9 #include "moments.hpp"
10 
12  int gid; // Corresponding XGC species
13  int grid_index; // Which grid this species is on
14  bool is_electron; // Whether species is electron
15  double T_avg;
16  double mass;
17  double mass_au;
18  double charge_eu;
20  double conv_factor;
21  double numeric_vth2; //[col_fm_core_init<col_fm_core]
22  double numeric_t; //[col_fm_core_init<col_fm_core]
23  double dens;
24  double ens;
25  double mom; // [col_fm_core_init<col_fm_core]
26 
27  //> Computes the base collision frequency -- Eq. (6) in Hager et al.
28  KOKKOS_INLINE_FUNCTION double lambda_gamma_pair(const CollisionSpeciesScalars& sp_b ) const {// C(a->b)
29  double lambda;
30  // NRL Plasma Formulary 2011 version, p.34
31  //(c) Mixed ion-ion collisions (here, same species ion-ion collisions)
32  if(is_electron || sp_b.is_electron){
33  //if electron is invovled
34  if(is_electron && sp_b.is_electron){
35  //(a) Thermal electron-electron collisions
36  // Note that pointers for colliding and target are same.
37  lambda = 2.35e1 - log( sqrt(den_lambda_gamma*1e-6)*pow(T_avg,-1.25))
38  - sqrt(std::abs(1e-5 + pow(log(T_avg)-2,2)/16));
39  } else {
40  double ti_ev,massi,densi,chargei,te_ev,masse,dense;
41  if(is_electron){
42  ti_ev = sp_b.T_avg;
43  massi = sp_b.mass_au;
44  densi = sp_b.den_lambda_gamma;
45  chargei = sp_b.charge_eu;
46  te_ev = T_avg;
47  masse = mass_au;
48  dense = den_lambda_gamma;
49  } else {
50  te_ev = sp_b.T_avg;
51  masse = sp_b.mass_au;
52  dense = sp_b.den_lambda_gamma;
53  ti_ev = T_avg;
54  massi = mass_au;
55  densi = den_lambda_gamma;
56  chargei = charge_eu;
57  }
58  //(b) Electron-ion (and ion-electron) collisions
59  if(ti_ev*masse/massi<te_ev){
60  if(te_ev<1e1*(chargei*chargei)){
61  lambda = 2.3e1 - log(chargei*sqrt(dense*1e-6/(te_ev*te_ev*te_ev)));
62  } else {
63  lambda = 2.4e1 - log(sqrt(dense*1e-6)/te_ev);
64  }
65  } else {
66  lambda = 3e1-log((chargei*chargei)*sqrt(densi*1e-6/(ti_ev*ti_ev*ti_ev))/massi);
67  }
68  }
69  } else {
70  lambda = 2.3e1 - log(charge_eu*sp_b.charge_eu*(mass_au+sp_b.mass_au)
71  /(mass_au*sp_b.T_avg+sp_b.mass_au*T_avg)
73  +sp_b.den_lambda_gamma*1e-6*(sp_b.charge_eu*sp_b.charge_eu)/sp_b.T_avg));
74  }
75  double tmp = lambda*(charge_eu*charge_eu)*(sp_b.charge_eu*sp_b.charge_eu)*pow(UNIT_CHARGE,4)/(pow(8.8542e-12,2)*8*PI);
76  return tmp/mass; // colliding : a, target : b
77  }
78 };
79 
80 template<class Device>
82  public:
84 
85  Kokkos::DualView<CollisionSpeciesScalars**,Device> s;
86 
87  Kokkos::DualView<double****,Kokkos::LayoutRight,Device> pdf_n; //[get_local_f<setup_col_sp] Given initial PDF
88  Kokkos::DualView<double****,Kokkos::LayoutRight,Device> pdf_np1; //[get_local_f<setup_col_sp] To be used as PDF at next Picard step. BUT has a role of "df" at end
89 
91 
92  CollisionSpecies<Device>(const VGridDistribution<HostType>& f_in, std::vector<Species<Device>> &all_species, int mb_n_nodes)
93  : f(f_in),
94  s("s", f.n_species(), mb_n_nodes),
95  pdf_n("pdf_n", mb_n_nodes, f.n_species(), f.n_vz(), f.n_vr()),
96  pdf_np1("pdf_np1", mb_n_nodes, f.n_species(), f.n_vz(), f.n_vr())
97  {
98  int col_species_cnt=0;
99  for (int isp = 0; isp<all_species.size(); isp++){
100  if (!all_species[isp].is_adiabatic){
101  for (int mesh_ind = 0; mesh_ind<mb_n_nodes; mesh_ind++){
102  s.h_view(col_species_cnt, mesh_ind).gid = isp;
103  s.h_view(col_species_cnt, mesh_ind).is_electron = all_species[isp].is_electron;
104  s.h_view(col_species_cnt, mesh_ind).grid_index = all_species[isp].collision_grid_index;
105  }
106  col_species_cnt++;
107  }
108  }
109  }
110 
111  //> Convert the distribution function to the units used
112  // internally by the collision operator
113  // The distribution function set up with f0_module is in normalized velocity space
114  // whereas the collision operator uses SI units.
115  void setup_one(int isp, int mesh_ind, const Species<Device> &species, const Moments& moments, int local_node_ind) {
116  s.h_view(isp,mesh_ind).mass = species.mass;
117  s.h_view(isp,mesh_ind).mass_au = s.h_view(isp,mesh_ind).mass/PROTON_MASS;
118  s.h_view(isp,mesh_ind).charge_eu = species.charge/UNIT_CHARGE;
119 
120  double conv_inv_mass=TWOPI*UNIT_CHARGE/s.h_view(isp,mesh_ind).mass;
121  s.h_view(isp,mesh_ind).den_lambda_gamma = moments.den_local(local_node_ind, species.nonadiabatic_idx);
122  s.h_view(isp,mesh_ind).T_avg = moments.temp_local(local_node_ind, species.nonadiabatic_idx);
123  double conv_fac_temp = species.f0.fg_temp_ev_h(local_node_ind);
124  s.h_view(isp,mesh_ind).conv_factor = 1.0/sqrt(conv_fac_temp*(conv_inv_mass*conv_inv_mass*conv_inv_mass));
125 
126  // Local f with basic correction:
127  // Simply cut off all negative values
128  for (int imu = 0; imu < pdf_n.h_view.extent(3); imu++){
129  for (int ivp = 0; ivp<pdf_n.h_view.extent(2); ivp++){
130  pdf_n.h_view(mesh_ind,isp,ivp,imu) = std::max(f(species.nonadiabatic_idx, imu,local_node_ind,ivp),0.0);
131  }
132  }
133 
134  // Convert to SI
135  for (int imu = 0; imu < pdf_n.h_view.extent(3); imu++){
136  double vnorm_to_SI = s.h_view(isp,mesh_ind).conv_factor/f.get_smu_n(imu);
137  for (int ivp = 0; ivp<pdf_n.h_view.extent(2); ivp++){
138  pdf_n.h_view(mesh_ind,isp,ivp,imu) *= vnorm_to_SI;
139  }
140  }
141 
142  // pdf_np1 = pdf_n - replace w deep copy - ALS
143  for (int iz = 0; iz < pdf_np1.h_view.extent(2); iz++){
144  for (int ir = 0; ir < pdf_np1.h_view.extent(3); ir++){
145  pdf_np1.h_view(mesh_ind,isp,iz,ir) = pdf_n.h_view(mesh_ind,isp,iz,ir);
146  }
147  }
148  }
149 
150  void setup_all(const std::vector<Species<Device>>& all_species, const Moments& moments, Kokkos::View<int**,HostType>& mesh_nodes, int mesh_batch_ind){
151 
152  for (int mesh_ind = 0; mesh_ind<s.extent(1); mesh_ind++){
153  for (int isp = 0; isp < s.extent(0); isp++){
154  setup_one(isp, mesh_ind, all_species[s.h_view(isp,mesh_ind).gid], moments, mesh_nodes(mesh_batch_ind,mesh_ind));
155  }
156  }
157 
158 #ifdef USE_GPU
159  // Copy dual views to device
160  Kokkos::deep_copy(s.d_view, s.h_view);
161  Kokkos::deep_copy(pdf_n.d_view, pdf_n.h_view);
162 #endif
163  }
164 
165  void lambda_gamma_pair(const View<double*,CLayout,HostType>& dt, Kokkos::View<double***,Device>& gammac){
166  // Done on host, so create a mirror
167  //typename Kokkos::View<double***,Device>::HostMirror
168  auto gammac_h = Kokkos::create_mirror_view(gammac);
169 
170  for (int mesh_ind = 0; mesh_ind<s.extent(1); mesh_ind++){
171  // gamma
172  for (int isp = 0; isp < s.extent(0); isp++){
173  for (int jsp = 0; jsp < s.extent(0); jsp++){
174  gammac_h(mesh_ind,jsp,isp) = dt(mesh_ind)*s.h_view(isp,mesh_ind).lambda_gamma_pair(s.h_view(jsp,mesh_ind));
175  }
176  }
177  }
178  Kokkos::deep_copy(gammac, gammac_h);
179  }
180 
181  inline int n() const{
182  return f.n_species();
183  }
184 };
185 
186 #endif
KOKKOS_INLINE_FUNCTION int n_vr() const
Definition: vgrid_distribution.hpp:155
constexpr double PROTON_MASS
Definition: constants.hpp:7
Distribution< Device > f0
Species distribution in velocity space on local mesh nodes.
Definition: species.hpp:127
View< double **, CLayout, HostType > temp_local
temperature, , local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:17
Kokkos::DualView< double ****, Kokkos::LayoutRight, Device > pdf_np1
Definition: col_species.hpp:88
double mass_au
Definition: col_species.hpp:17
double get_smu_n(int imu) const
Definition: vgrid_distribution.hpp:172
double den_lambda_gamma
Definition: col_species.hpp:19
bool is_electron
Definition: col_species.hpp:14
Definition: moments.hpp:13
int gid
Definition: col_species.hpp:12
void setup_one(int isp, int mesh_ind, const Species< Device > &species, const Moments &moments, int local_node_ind)
Definition: col_species.hpp:115
double conv_factor
Definition: col_species.hpp:20
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays) ...
Definition: species.hpp:81
Kokkos::DualView< CollisionSpeciesScalars **, Device > s
Definition: col_species.hpp:85
Definition: col_species.hpp:11
KOKKOS_INLINE_FUNCTION int n_species() const
Definition: vgrid_distribution.hpp:151
double mass
Particle mass.
Definition: species.hpp:83
double numeric_t
Definition: col_species.hpp:22
KOKKOS_INLINE_FUNCTION int n_vz() const
Definition: vgrid_distribution.hpp:163
Kokkos::DualView< double ****, Kokkos::LayoutRight, Device > pdf_n
Definition: col_species.hpp:87
void lambda_gamma_pair(const View< double *, CLayout, HostType > &dt, Kokkos::View< double ***, Device > &gammac)
Definition: col_species.hpp:165
double mass
Definition: col_species.hpp:16
double charge
Particle charge.
Definition: species.hpp:84
View< double **, CLayout, HostType > den_local
density, local to compute node [Nsp,Nlocalnode]
Definition: moments.hpp:14
double charge_eu
Definition: col_species.hpp:18
void setup_all(const std::vector< Species< Device >> &all_species, const Moments &moments, Kokkos::View< int **, HostType > &mesh_nodes, int mesh_batch_ind)
Definition: col_species.hpp:150
constexpr double UNIT_CHARGE
Charge of an electron (C)
Definition: constants.hpp:4
double numeric_vth2
Definition: col_species.hpp:21
int grid_index
Definition: col_species.hpp:13
int n() const
Definition: col_species.hpp:181
double T_avg
Definition: col_species.hpp:15
double dens
Definition: col_species.hpp:23
double mom
Definition: col_species.hpp:25
const VGridDistribution< HostType > f
Definition: col_species.hpp:83
Definition: species.hpp:75
Definition: col_species.hpp:81
constexpr double PI
Definition: constants.hpp:8
constexpr double TWOPI
Definition: constants.hpp:9
double ens
Definition: col_species.hpp:24
KOKKOS_INLINE_FUNCTION double lambda_gamma_pair(const CollisionSpeciesScalars &sp_b) const
Definition: col_species.hpp:28