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