1 #ifndef MONTE_CARLO_COLLISIONS_HPP
2 #define MONTE_CARLO_COLLISIONS_HPP
34 vb.update(magnetic_field, species);
55 RandGen rand_gen = rand_pool.get_state();
58 for (
int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
62 double rho_m = part.
ph.
rho[i_simd];
66 double mu = part.
ct.
mu[i_simd];
70 ekin = max(ekmin,ekin);
74 double vprt = sqrt(2.0*ekin/species.
mass);
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);
82 part.
ph.
rho[i_simd] = rho_m;
83 part.
ct.
mu[i_simd] = mu;
87 rand_pool.free_state(rand_gen);
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);
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;
118 double del_pitch = random_sign*sqrt(dum*colb*dt*0.5);
119 pitch = pitch*(1.0 - colb*dt*0.5) + del_pitch;
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);
129 DEVICE_PRINTF(
"\nERROR: ekin is nan: %1.15e; colbs=%1.15e, fac0b=%1.15e, esig=%1.15e \n", ekin, colbs, fac0b, esig);
133 pitch = min(1.0, pitch);
134 pitch = max(-1.0, pitch);
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{
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;
166 double cnst = 2.41e5*charge*charge/(vprt*vprt*vprt*mass*mass);
172 double dn = dn_b/1.0e6;
176 double ee_ev = en_a*
J_2_EV;
179 double dumb = vprt/vt;
180 double dd = dumb*dumb;
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);
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;
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);
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));
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);
219 up=(ftype>=1)? up*r : up;
224 up *= bphi/sqrt(br*br+bz*bz+bphi*bphi);
236 accel = nlr.
get<
bool>(
"col_accel",
false);
243 double psi_range = magnetic_field.
outpsi - magnetic_field.
inpsi;
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");
264 if(
vb.period<
period)
exit_XGC(
"\nError: must set col_vb_period >= col_period\n");
271 TIMER(
"copy_ptl_to_device",
275 TIMER(
"COL_SNAPSHOT",
276 vb.update(magnetic_field, species) );
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 ¶m, const T default_val, int val_ind=0)
Definition: NamelistReader.hpp:372
bool en_col_on
Whether to use energy collisions.
Definition: monte_carlo_collisions.hpp:19
double c2_2m
c2/2m
Definition: species.hpp:85
#define DEVICE_PRINTF(...)
Definition: space_settings.hpp:85
double c_m
c/m
Definition: species.hpp:84
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:123
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:584
Definition: NamelistReader.hpp:193
Definition: magnetic_field.hpp:12
void for_all_nonadiabatic_species(F func, DevicePtlOpt device_ptl_opt=UseDevicePtl)
Definition: plasma.hpp:163
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:79
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:125
double charge_eu
Particle charge in eu.
Definition: species.hpp:83
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:81
#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:339
void use_namelist(const string &namelist)
Definition: NamelistReader.hpp:354
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:124
SimdConstants ct
Definition: particles.hpp:60
Definition: plasma.hpp:15
int period
Definition: monte_carlo_collisions.hpp:27
Definition: species.hpp:73
double accel_pout2
Definition: monte_carlo_collisions.hpp:16
Eq::Profile< Device > eq_temp
Definition: species.hpp:122
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