1 #ifndef MONTE_CARLO_COLLISIONS_HPP
2 #define MONTE_CARLO_COLLISIONS_HPP
36 vb.update(magnetic_field, species);
57 RandGen rand_gen = rand_pool.get_state();
60 for (
int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
64 double rho_m = part.
ph.
rho[i_simd];
68 double mu = part.
ct.
mu[i_simd];
72 ekin = max(ekmin,ekin);
76 double vprt = sqrt(2.0*ekin/species.
mass);
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);
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;
91 rand_pool.free_state(rand_gen);
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);
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;
122 double del_pitch = random_sign*sqrt(dum*colb*dt*0.5);
123 pitch = pitch*(1.0 - colb*dt*0.5) + del_pitch;
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);
133 DEVICE_PRINTF(
"\nERROR: ekin is nan: %1.15e; colbs=%1.15e, fac0b=%1.15e, esig=%1.15e \n", ekin, colbs, fac0b, esig);
137 pitch = min(1.0, pitch);
138 pitch = max(-1.0, pitch);
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{
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;
170 double cnst = 2.41e5*charge*charge/(vprt*vprt*vprt*mass*mass);
176 double dn = dn_b/1.0e6;
180 double ee_ev = en_a*
J_2_EV;
183 double dumb = vprt/vt;
184 double dd = dumb*dumb;
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);
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;
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);
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));
218 den = species.
eq_den.
value(magnetic_field, psi, r, z);
219 temp = species.
eq_temp.
value(magnetic_field, psi, r, z);
223 up=(ftype>=1)? up*r : up;
228 up *= bphi/sqrt(br*br+bz*bz+bphi*bphi);
240 accel = nlr.
get<
bool>(
"col_accel",
false);
247 double psi_range = magnetic_field.
outpsi - magnetic_field.
inpsi;
252 pin = nlr.
get<
double>(
"col_pin" , magnetic_field.
inpsi );
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");
275 if(
vb.period<
period)
exit_XGC(
"\nError: must set col_vb_period >= col_period\n");
282 TIMER(
"copy_ptl_to_device",
286 TIMER(
"COL_SNAPSHOT",
287 vb.update(magnetic_field, species) );
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: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:113
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:204
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: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 ¶m, 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 ¶m)
Definition: NamelistReader.hpp:363
Eq::Profile< Device > eq_den
Definition: species.hpp:134
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:600
Definition: NamelistReader.hpp:193
Definition: magnetic_field.hpp:12
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: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:142
KOKKOS_INLINE_FUNCTION void bmag_interpol(const SimdVector &v, Simd< double > &bmag) const
Definition: magnetic_field.tpp:239
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 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:136
double charge_eu
Particle charge in eu.
Definition: species.hpp:85
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector2D &x, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:57
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:26
void copy_particles_to_device_if_not_resident()
Definition: species.hpp:362
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:129
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:135
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:133
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
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:79
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