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