XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
monte_carlo_collisions.hpp
Go to the documentation of this file.
1 #ifndef MONTE_CARLO_COLLISIONS_HPP
2 #define MONTE_CARLO_COLLISIONS_HPP
3 
4 #include "basic_physics.hpp"
5 #include "varying_background.hpp"
6 #include "magnetic_field.hpp"
7 #include "plasma.hpp"
8 
9 template<class Device>
11  bool accel;
12  int accel_n;
13  double accel_pin1;
14  double accel_pout1;
15  double accel_pin2;
16  double accel_pout2;
17  double accel_factor1;
18  double accel_factor2;
19  bool en_col_on;
21  double pin;
22  double pout;
23 
26 
27  public:
28 
29  int period;
30 
31  inline bool do_snapshot(int istep) const{
32  return (use_varying_bg && (istep==1 || (istep%vb.period==0)));
33  }
34 
36  vb.update(magnetic_field, species);
37  }
38 
39  KOKKOS_INLINE_FUNCTION void collision_c(const MagneticField<DeviceType>& magnetic_field, const Species<DeviceType>& species, const Species<DeviceType>& species_b, const pool_type& rand_pool, const double dt, const bool do_pitch_angle_scattering, const double ekmin, const int idx) const{
40  // Get index info
41  AoSoAIndices<DeviceType> inds(idx);
42 
43  // Copy from AoSoA to SoA
44  SimdParticles part;
45  AoSoA_2_simd(part, species.ptl(), inds.a, inds.s);
46 
47  Simd<double> psi;
48  magnetic_field.get_psi(part.ph.v(),psi);
49 
50  Simd<double> Bmag;
51  magnetic_field.bmag_interpol(part.ph.v(), Bmag);
52 
53  Simd<double> theta;
54  magnetic_field.equil.get_theta(part.ph.x(), theta);
55 
56  // Sample random numbers for each accumulator; double is totally overkill since we just use it for a sign
57  RandGen rand_gen = rand_pool.get_state();
58 
60  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
61  double dnb, t_ev, up;
62  background_sp_profile(magnetic_field, theta[i_simd],part.ph.r[i_simd],part.ph.z[i_simd],psi[i_simd],species_b, dnb,t_ev,up);
63 
64  double rho_m = part.ph.rho[i_simd];
65  if(moving_frame) rho_m -= up/(species_b.c_m*Bmag[i_simd]);
66 
67  //test ptl collision
68  double mu = part.ct.mu[i_simd];
69 
70  double pitch, ekin;
71  rho_mu_to_en_pitch(Bmag[i_simd], species.c_m, species.c2_2m, species.mass, rho_m, mu, ekin, pitch);
72  ekin = max(ekmin,ekin); // for preventing small energy
73 
74  double accel_factor = (accel ? get_accel_factor(psi[i_simd]) : 1.0);
75 
76  double vprt = sqrt(2.0*ekin/species.mass); // MKS - m/s
77  scatter_one(rand_gen, vprt, species.mass/PROTON_MASS,species.charge_eu,dnb,t_ev,species_b.mass/PROTON_MASS,species_b.charge_eu,accel_factor,dt,ekmin,do_pitch_angle_scattering, ekin, pitch);
78 
79  en_pitch_to_rho_mu(Bmag[i_simd], species.c_m, species.mass, ekin, pitch, rho_m, mu);
80 
81  if(moving_frame) rho_m += up/(species_b.c_m*Bmag[i_simd]);
82 
83  // Update particle
84  if (psi[i_simd] >= pin && psi[i_simd] <= pout){
85  part.ph.rho[i_simd] = rho_m;
86  part.ct.mu[i_simd] = mu;
87  }
88  }
89 
90  // Give the state back, which will allow another thread to acquire it
91  rand_pool.free_state(rand_gen);
92 
93  simd_2_AoSoA(species.ptl(), part, inds.a, inds.s);
94  }
95 
96  private:
97 
98  KOKKOS_INLINE_FUNCTION double get_accel_factor(double psi) const{
99  double factor = 1.0;
100  if(accel_n>=1){
101  if(accel_pin1 < psi && psi < accel_pout1 ){
102  factor *= accel_factor1;
103  }
104  }
105  if(accel_n>=2){
106  if(accel_pin2 < psi && psi < accel_pout2 ){
107  factor *= accel_factor2;
108  }
109  }
110  return factor;
111  }
112 
113  KOKKOS_INLINE_FUNCTION void scatter_one(RandGen& rand_gen, double vprt, double massa_au, double chargea_eu, double denb, double tempb_ev, double massb_au, double chargeb_eu, double accel, double dt, double ekmin, bool do_pitch_angle_scattering, double& ekin, double& pitch) const {
114  double colb, colbs, fac0b;
115  find_freq(ekin,vprt, massa_au,chargea_eu,denb,tempb_ev,massb_au,chargeb_eu,accel, colb,colbs,fac0b);
116 
117  // pitch angle scattering
118  if(do_pitch_angle_scattering){
119  double random_sign = (rand_gen.drand() >= 0.5 ? 1.0 : -1.0);
120  double dum = 1.0 - pitch*pitch;
121  dum = max(0.0, dum);
122  double del_pitch = random_sign*sqrt(dum*colb*dt*0.5);
123  pitch = pitch*(1.0 - colb*dt*0.5) + del_pitch;
124  }
125 
126  // energy collision
127  if(en_col_on){
128  double random_sign = (rand_gen.drand() >= 0.5 ? 1.0 : -1.0);
129  double esig = random_sign*sqrt(2*ekin*tempb_ev*EV_2_J*colbs*dt);
130  ekin = ekin - colbs*dt*fac0b + esig;
131  ekin = max(ekin, ekmin);
132  if(isnan(ekin)){
133  DEVICE_PRINTF("\nERROR: ekin is nan: %1.15e; colbs=%1.15e, fac0b=%1.15e, esig=%1.15e \n", ekin, colbs, fac0b, esig);
134  }
135  }
136 
137  pitch = min(1.0, pitch);
138  pitch = max(-1.0, pitch);
139  }
140 
141  // Warning : This routine uses cgs unit partially
142  KOKKOS_INLINE_FUNCTION void find_freq(double en_a, double vprt, double mass, double charge, double dn_b, double en_b_ev, double mass_b, double charge_b, double accel, double& freq_scat, double& freq_slow, double& freq_fac0) const{
143  // calculate collision frequencies and convert into MKS unit
144  // approx to psi function
145  // pitch angle scattering, small value of col! (alpha)
146  // Boozer and Kuo-Petravic Phys. Fluids 24, 851 (1981)
147  // profiles-density (cm-3) dnb (background), dni (impurity), temp (kev)
148  // psi(x) = (2/sqrt(pi)Int[dt*sqrt(t)*exp(-t)]
149  // psi(x) = Phi(v/vth) = Phi(sqrt(x)) in the paper
150  // psi(x) = 1 - 1/(1 + p), approximatly
151  // p = x**1.5[a0 + a1*x +a2*x**2 +a3*x**3 +a4*x**4]
152  // R.White,M.Redi Jan (2000) error in psi(x) < 1.D-3, asymptotically correct
153  // relative error dpsi/psi < 1.D-2
154  // en_a(J), mass(kg), charge(C), dn_b (m^-3), en_b_ev (eV), mass_b (Atomic Unit), charge_b (Electron Charge unit)
155 
156  constexpr double ap0 = .75225;
157  constexpr double ap1 = -.238;
158  constexpr double ap2 = .311;
159  constexpr double ap3 = -.0956;
160  constexpr double ap4 = .0156;
161 
162  // vt = dsqrt(2d0*(en_b_ev/1.d3*sml_en_order/sml_en_order_kev)*mass/mass_b)
163  double vt = sqrt(2.*(en_b_ev*EV_2_J)/(mass_b*PROTON_MASS)); // MKS - m/s
164 
165  // norm_r_cgs = nc_norm_r * 100.d0 !!! conversion MKS to cgs
166  // cnst = 2.41D11*(charge/mass)**2/(norm_r_cgs/nc_norm_t*vprt)**3 !100 is conversion MKS to CGS
167  // cnst = cnst*col_accel
168 
169  // 2.41e5 = 2.41e11/1.0e6 //1e2^3 is to convert vprt from MKS to CGS
170  double cnst = 2.41e5*charge*charge/(vprt*vprt*vprt*mass*mass);
171  // cnst = mu_b / x^3 /Clog /den / some order 1 const
172  // mu_b=Braginskii collision frequency
173  // x = v/vth, Clog = Coulumb Logarithm
174  // 2.4D11 = e(cgs)^4 / m(cgs)^2 * 4*pi
175  cnst *= accel;
176  double dn = dn_b/1.0e6; // dn - CGS cm^-3
177 
178 
179 
180  double ee_ev = en_a*J_2_EV; // test ptl's kinetic energy in [eV]
181 
182  // calculate psi(x), f, g, gp
183  double dumb = vprt/vt; // MKS/MKS , v/v_th
184  double dd = dumb*dumb; // dd = x, dumb = sqrt(x)
185  double d3 = dumb*dd;
186  double dum = d3*(ap0 + ap1*dd + ap2*dd*dd + ap3*dd*dd*dd + ap4*dd*dd*dd*dd);
187  double ap_psi = 1.0 - 1.0/(1.0 + dum); //psi(x) = Phi(v/vth)
188 
189  double f = (2.0 - 1.0/dd)*ap_psi+2.257*dumb*exp(-dd);
190  double g = 2.0*mass*ap_psi/mass_b;
191  double gp = 2.257*exp(-dd)*(mass/mass_b)*dumb;
192 
193  // find Coulomb logarithm
194  // bmax = 7.4d-3*dsqrt(en_b_ev/1000d0*1.d13/(charge**2*dn_b))
195  // bmin1 = 1.4d-10*charge*charge_b/(ee_ev + en_b_ev)*1000d0
196  // massab = mass_au*mass_b/(mass_au+mass_b)
197  // bmin2 = 6.2d-10/(dsqrt(massab*1836.d0*en_b_ev/1000d0))
198 
199 
200  double bmax = 7.4e-3*sqrt(en_b_ev*1.e10/(charge*charge*dn));
201  double bmin1 = 1.4e-7*abs(charge*charge_b)/(ee_ev+en_b_ev);
202  double massab = mass*mass_b/(mass + mass_b);
203  double bmin2 = 6.2e-10/(sqrt(massab*1836.*en_b_ev/1.e3));
204  double bmin = max(bmin1,bmin2);
205  double clog = log(bmax/bmin);
206 
207  // collision frequencies - scattering, slowing down
208  double dnu_b = cnst*dn*charge_b*charge_b;
209  freq_scat = abs(clog*dnu_b*f);
210  freq_slow = abs(clog*dnu_b*g);
211  freq_fac0 = en_a*(1.0 - mass_b*gp/(g*mass));
212  }
213 
214  KOKKOS_INLINE_FUNCTION void background_sp_profile(const MagneticField<DeviceType>& magnetic_field, double theta, double r, double z, double psi, const Species<DeviceType>& species, double& den, double& temp, double& up) const{
215  if(use_varying_bg && vb.is_in_range(magnetic_field, r, z, psi)){
216  vb.interp_bg_profile(psi,theta,species.nonadiabatic_idx, den, temp, up);
217  } else {
218  den = species.eq_den.value(magnetic_field, psi, r, z);
219  temp = species.eq_temp.value(magnetic_field, psi, r, z);
220  up = species.eq_flow.value(magnetic_field, psi, r, z);
221 
222  int ftype=species.eq_flow_type;
223  up=(ftype>=1)? up*r : up;
224  if(ftype>=2) { // ftype==2 might slow down the performance a litte.
225  double phi=0; // assign random angle assuming toroidal symmetry.
226  double br,bz,bphi;
227  magnetic_field.bvec_interpol(r,z,phi,br,bz,bphi);
228  up *= bphi/sqrt(br*br+bz*bz+bphi*bphi);
229  }
230  }
231  }
232 
233  public:
234 
236 
238  nlr.use_namelist("col_param");
239  period = nlr.get<int>("col_period", 3);
240  accel = nlr.get<bool>("col_accel", false);
241 
242  accel_n = nlr.get<int>("col_accel_n",2);
243  accel_factor1 = nlr.get<double>("col_accel_factor1",10.0);
244  accel_factor2 = nlr.get<double>("col_accel_factor2",10.0);
245 
246  accel_pin1 = nlr.get<double>("col_accel_pin1", magnetic_field.inpsi);
247  double psi_range = magnetic_field.outpsi - magnetic_field.inpsi;
248  accel_pout1 = nlr.get<double>("col_accel_pout1", magnetic_field.inpsi + 0.1*psi_range);
249  accel_pin2 = nlr.get<double>("col_accel_pin2", magnetic_field.outpsi - 0.1*psi_range);
250  accel_pout2 = nlr.get<double>("col_accel_pout2", magnetic_field.outpsi);
251 
252  pin = nlr.get<double>("col_pin" , magnetic_field.inpsi );
253  pout = nlr.get<double>("col_pout", magnetic_field.outpsi);
254 
255  // Normalize
256  accel_pin1 *= magnetic_field.psi_norm();
257  accel_pout1 *= magnetic_field.psi_norm();
258  accel_pin2 *= magnetic_field.psi_norm();
259  accel_pout2 *= magnetic_field.psi_norm();
260  // Normalize by xpt
261  // Default values are already normalized, so only normalize if not default
262  if(nlr.present("col_pin")) pin *= magnetic_field.psi_norm();
263  if(nlr.present("col_pout")) pin *= magnetic_field.psi_norm();
264 
265  en_col_on = (nlr.get<int>("col_en_col_on", 1))==1;
266  moving_frame = nlr.get<bool>("col_moving_frame", true);
267 
268  use_varying_bg = (nlr.get<int>("col_varying_bg", 0))==1;
269  if(use_varying_bg){
270  exit_XGC("\nThe option col_varying_bg is not currently operational. It will work only in a full-f simulation. Delta-f does not need a varying background. But for total-f, this routine would have to be extended/modified.\n");
271 
272  // Initialize varying background
273  vb = VaryingBackground<Device>(nlr, pol_decomp, magnetic_field, grid, n_nonadiabatic_species);
274 
275  if(vb.period<period) exit_XGC("\nError: must set col_vb_period >= col_period\n");
276  }
277  }
278 
281  // Copy particles to the device
282  TIMER("copy_ptl_to_device",
284 
285  // A snapshot is created for ions and electrons every col_vb_period calls to collision
286  TIMER("COL_SNAPSHOT",
287  vb.update(magnetic_field, species) );
288  });
289  }
290 };
291 
293 
294 #endif
Kokkos::Random_XorShift64_Pool< DeviceType > pool_type
Definition: space_settings.hpp:71
double inpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:43
KOKKOS_INLINE_FUNCTION void scatter_one(RandGen &rand_gen, double vprt, double massa_au, double chargea_eu, double denb, double tempb_ev, double massb_au, double chargeb_eu, double accel, double dt, double ekmin, bool do_pitch_angle_scattering, double &ekin, double &pitch) const
Definition: monte_carlo_collisions.hpp:113
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:321
constexpr double PROTON_MASS
Definition: constants.hpp:7
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector &v, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:105
KOKKOS_INLINE_FUNCTION void background_sp_profile(const MagneticField< DeviceType > &magnetic_field, double theta, double r, double z, double psi, const Species< DeviceType > &species, double &den, double &temp, double &up) const
Definition: monte_carlo_collisions.hpp:214
KOKKOS_INLINE_FUNCTION void simd_2_AoSoA(VecParticles *part, const SimdParticles &part_one, int a_vec, int s_vec)
Definition: particles.tpp:54
T get(const string &param, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:373
bool en_col_on
Whether to use energy collisions.
Definition: monte_carlo_collisions.hpp:19
double c2_2m
c2/2m
Definition: species.hpp:87
#define DEVICE_PRINTF(...)
Definition: space_settings.hpp:85
double c_m
c/m
Definition: species.hpp:86
constexpr double EV_2_J
Conversion rate ev to J.
Definition: constants.hpp:5
double accel_factor2
Definition: monte_carlo_collisions.hpp:18
bool moving_frame
Definition: monte_carlo_collisions.hpp:20
bool present(const string &param)
Definition: NamelistReader.hpp:363
Eq::Profile< Device > eq_den
Definition: species.hpp:132
subroutine plasma(grid, itr, p, dene_out, deni_out, Te_out, Ti_out, Vparai_out)
Calculate the plasma density, temperature, and parallel velocity for a point in triangle itr using pl...
Definition: neutral_totalf.F90:1224
KOKKOS_INLINE_FUNCTION VecParticles * ptl() const
Definition: species.hpp:597
Definition: NamelistReader.hpp:193
Definition: magnetic_field.hpp:13
void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:132
int a
The index in the inner array of the AoSoA.
Definition: particles.hpp:150
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:45
KOKKOS_INLINE_FUNCTION void find_freq(double en_a, double vprt, double mass, double charge, double dn_b, double en_b_ev, double mass_b, double charge_b, double accel, double &freq_scat, double &freq_slow, double &freq_fac0) const
Definition: monte_carlo_collisions.hpp:142
KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector &v, Simd< double > &bmag) const
Definition: magnetic_field.tpp:372
double pout
Outer boundary for collisions in norm. pol. flux.
Definition: monte_carlo_collisions.hpp:22
int nonadiabatic_idx
Index of species skipping adiabatic species (for compatibility with fortran arrays) ...
Definition: species.hpp:81
KOKKOS_INLINE_FUNCTION void en_pitch_to_rho_mu(double B, double c_m, double mass, double en, double pitch, double &rho, double &mu)
Definition: basic_physics.hpp:103
KOKKOS_INLINE_FUNCTION double psi_norm() const
Definition: magnetic_field.tpp:62
KOKKOS_INLINE_FUNCTION void AoSoA_2_simd(SimdParticles &part_one, const VecParticles *part, int a_vec, int s_vec)
Definition: particles.tpp:80
Simd< double > rho
Definition: particles.hpp:21
double accel_factor1
Definition: monte_carlo_collisions.hpp:17
double accel_pin2
Definition: monte_carlo_collisions.hpp:15
int eq_flow_type
Definition: species.hpp:134
double charge_eu
Particle charge in eu.
Definition: species.hpp:85
double accel_pout1
Definition: monte_carlo_collisions.hpp:14
constexpr double J_2_EV
Conversion rate J to ev.
Definition: constants.hpp:6
double mass
Particle mass.
Definition: species.hpp:83
#define TIMER(N, F)
Definition: timer_macro.hpp:24
double outpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:44
void copy_particles_to_device_if_not_resident()
Definition: species.hpp:359
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:355
Simd< double > r
Definition: particles.hpp:18
bool use_varying_bg
Definition: monte_carlo_collisions.hpp:24
KOKKOS_INLINE_FUNCTION SimdVector2D & x()
Definition: particles.hpp:27
KOKKOS_INLINE_FUNCTION double get_accel_factor(double psi) const
Definition: monte_carlo_collisions.hpp:98
KOKKOS_INLINE_FUNCTION SimdVector & v()
Definition: particles.hpp:39
KOKKOS_INLINE_FUNCTION void rho_mu_to_en_pitch(double B, double c_m, double c2_2m, double mass, double rho, double mu, double &en, double &pitch)
Definition: basic_physics.hpp:95
SimdPhase ph
Definition: particles.hpp:59
pool_type::generator_type RandGen
Definition: space_settings.hpp:72
Definition: particles.hpp:58
void update_vb(const MagneticField< DeviceType > &magnetic_field, Species< DeviceType > &species)
Definition: monte_carlo_collisions.hpp:35
void monte_carlo_collisions(int istep, const MagneticField< DeviceType > &magnetic_field, Plasma &plasma, MonteCarloCollider< DeviceType > &col_mc, double dt_in)
Definition: monte_carlo_collisions.cpp:22
int s
The index in the outer array of the AoSoA.
Definition: particles.hpp:149
#define FORCE_SIMD
Definition: simd.hpp:9
Simd< double > z
Definition: particles.hpp:19
int accel_n
Definition: monte_carlo_collisions.hpp:12
KOKKOS_INLINE_FUNCTION void get_theta(const SimdVector2D &x, Simd< double > &theta) const
Definition: equil.tpp:142
bool do_snapshot(int istep) const
Definition: monte_carlo_collisions.hpp:31
void exit_XGC(std::string msg)
Definition: globals.hpp:37
Definition: magnetic_field.F90:1
KOKKOS_INLINE_FUNCTION void collision_c(const MagneticField< DeviceType > &magnetic_field, const Species< DeviceType > &species, const Species< DeviceType > &species_b, const pool_type &rand_pool, const double dt, const bool do_pitch_angle_scattering, const double ekmin, const int idx) const
Definition: monte_carlo_collisions.hpp:39
Eq::Profile< Device > eq_flow
Definition: species.hpp:133
SimdConstants ct
Definition: particles.hpp:60
KOKKOS_INLINE_FUNCTION double value(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:319
Definition: plasma.hpp:14
double pin
Inner boundary for collisions in norm. pol. flux.
Definition: monte_carlo_collisions.hpp:21
int period
Definition: monte_carlo_collisions.hpp:29
Definition: species.hpp:75
double accel_pout2
Definition: monte_carlo_collisions.hpp:16
Eq::Profile< Device > eq_temp
Definition: species.hpp:131
double accel_pin1
Definition: monte_carlo_collisions.hpp:13
Simd< double > mu
Definition: particles.hpp:52
Definition: varying_background.hpp:10
Definition: particles.hpp:148
void update_vb_all_species(const MagneticField< DeviceType > &magnetic_field, Plasma &plasma)
Definition: monte_carlo_collisions.hpp:279
MonteCarloCollider(NLReader::NamelistReader &nlr, const DomainDecomposition< DeviceType > &pol_decomp, const MagneticField< DeviceType > &magnetic_field, const Grid< DeviceType > &grid, int n_nonadiabatic_species)
Definition: monte_carlo_collisions.hpp:237
Definition: monte_carlo_collisions.hpp:10
bool accel
Artificial acceleration of collisions.
Definition: monte_carlo_collisions.hpp:11
VaryingBackground< Device > vb
Definition: monte_carlo_collisions.hpp:25
MonteCarloCollider()
Definition: monte_carlo_collisions.hpp:235