monte_carlo_collisions.hpp

Go to the documentation of this file.

#ifndef MONTE_CARLO_COLLISIONS_HPP
#define MONTE_CARLO_COLLISIONS_HPP

#include "basic_physics.hpp"
#include "varying_background.hpp"
#include "magnetic_field.hpp"
#include "plasma.hpp"

template<class Device>
class MonteCarloCollider{
    bool accel; 
    int accel_n;
    double accel_pin1;
    double accel_pout1;
    double accel_pin2;
    double accel_pout2;
    double accel_factor1;
    double accel_factor2;
    bool en_col_on; 
    bool moving_frame;

    bool use_varying_bg;
    VaryingBackground<Device> vb;

    public:

    int period;
   
    inline bool do_snapshot(int istep) const{
        return (use_varying_bg && (istep==1 || (istep%vb.period==0)));
    }

    void update_vb(const MagneticField<DeviceType>& magnetic_field, Species<DeviceType>& species){
        vb.update(magnetic_field, species);
    }

    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{
        // Get index info
        AoSoAIndices<DeviceType> inds(idx);

        // Copy from AoSoA to SoA
        SimdParticles part;
        AoSoA_2_simd(part, species.ptl(), inds.a, inds.s);

        Simd<double> psi;
        magnetic_field.get_psi(part.ph.x(),psi);

        Simd<double> Bmag;
        magnetic_field.bmag_interpol(part.ph.v(), Bmag);

        Simd<double> theta;
        magnetic_field.equil.get_theta(part.ph.x(), theta);

        // Sample random numbers for each accumulator; double is totally overkill since we just use it for a sign
        RandGen rand_gen = rand_pool.get_state();

        FORCE_SIMD
        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            double dnb, t_ev, up;
            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);

            double rho_m = part.ph.rho[i_simd];
            if(moving_frame) rho_m -= up/(species_b.c_m*Bmag[i_simd]);

            //test ptl collision
            double mu = part.ct.mu[i_simd];

            double pitch, ekin;
            rho_mu_to_en_pitch(Bmag[i_simd], species.c_m, species.c2_2m, species.mass, rho_m,  mu, ekin, pitch);
            ekin = max(ekmin,ekin); // for preventing small energy

            double accel_factor = (accel ? get_accel_factor(psi[i_simd]) : 1.0);

            double vprt = sqrt(2.0*ekin/species.mass); // MKS - m/s
            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);

            en_pitch_to_rho_mu(Bmag[i_simd], species.c2_2m, ekin, pitch, rho_m, mu);

            if(moving_frame) rho_m += up/(species_b.c_m*Bmag[i_simd]);

            // Update particle
            part.ph.rho[i_simd] = rho_m;
            part.ct.mu[i_simd] = mu;
        }

        // Give the state back, which will allow another thread to acquire it
        rand_pool.free_state(rand_gen);

        simd_2_AoSoA(species.ptl(), part, inds.a, inds.s);
    }

    private:

    KOKKOS_INLINE_FUNCTION double get_accel_factor(double psi) const{
        double factor = 1.0;
        if(accel_n>=1){
            if(accel_pin1 < psi && psi < accel_pout1 ){
                factor *= accel_factor1;
            }
        }
        if(accel_n>=2){
            if(accel_pin2 < psi && psi < accel_pout2 ){
                factor *= accel_factor2;
            }
        }
        return factor;
    }

    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 {
        double colb, colbs, fac0b;
        find_freq(ekin,vprt, massa_au,chargea_eu,denb,tempb_ev,massb_au,chargeb_eu,accel, colb,colbs,fac0b);

        // pitch angle scattering
        if(do_pitch_angle_scattering){
            double random_sign = (rand_gen.drand() >= 0.5 ? 1.0 : -1.0);
            double dum = 1.0 - pitch*pitch;
            dum = max(0.0, dum);
            double del_pitch = random_sign*sqrt(dum*colb*dt*0.5);
            pitch = pitch*(1.0 - colb*dt*0.5) + del_pitch;
        }

        // energy collision 
        if(en_col_on){
            double random_sign = (rand_gen.drand() >= 0.5 ? 1.0 : -1.0);
            double esig = random_sign*sqrt(2*ekin*tempb_ev*EV_2_J*colbs*dt);
            ekin = ekin - colbs*dt*fac0b + esig;
            ekin = max(ekin, ekmin);
            if(isnan(ekin)){
                DEVICE_PRINTF("\nERROR: ekin is nan: %1.15e; colbs=%1.15e, fac0b=%1.15e, esig=%1.15e \n", ekin, colbs, fac0b, esig);
            }
        }

        pitch = min(1.0, pitch);
        pitch = max(-1.0, pitch);
    }

    // Warning : This routine uses cgs unit partially
    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{
        // calculate collision frequencies and convert into MKS unit
        // approx to psi function
        // pitch angle scattering, small value of col!  (alpha)
        // Boozer and Kuo-Petravic Phys. Fluids 24, 851 (1981)
        // profiles-density (cm-3) dnb (background), dni (impurity), temp (kev)
        // psi(x) = (2/sqrt(pi)Int[dt*sqrt(t)*exp(-t)]
        // psi(x) = Phi(v/vth) = Phi(sqrt(x)) in the paper
        // psi(x) = 1 - 1/(1 + p), approximatly
        // p = x**1.5[a0 + a1*x +a2*x**2 +a3*x**3 +a4*x**4]
        // R.White,M.Redi Jan (2000) error in psi(x) < 1.D-3, asymptotically correct
        // relative error dpsi/psi < 1.D-2
        // en_a(J), mass(kg), charge(C), dn_b (m^-3), en_b_ev (eV), mass_b (Atomic Unit), charge_b (Electron Charge unit)

        constexpr double ap0 = .75225;
        constexpr double ap1 = -.238;
        constexpr double ap2 = .311;
        constexpr double ap3 = -.0956;
        constexpr double ap4 = .0156;

        //  vt = dsqrt(2d0*(en_b_ev/1.d3*sml_en_order/sml_en_order_kev)*mass/mass_b)
        double vt = sqrt(2.*(en_b_ev*EV_2_J)/(mass_b*PROTON_MASS)); // MKS - m/s

        //  norm_r_cgs = nc_norm_r * 100.d0 !!! conversion MKS to cgs
        //  cnst = 2.41D11*(charge/mass)**2/(norm_r_cgs/nc_norm_t*vprt)**3  !100 is conversion MKS to CGS
        //  cnst = cnst*col_accel

        //  2.41e5 = 2.41e11/1.0e6 //1e2^3 is to convert vprt from MKS to CGS
        double cnst = 2.41e5*charge*charge/(vprt*vprt*vprt*mass*mass);
        // cnst = mu_b / x^3 /Clog /den / some order 1 const
        // mu_b=Braginskii collision frequency
        // x = v/vth, Clog = Coulumb Logarithm
        // 2.4D11 = e(cgs)^4 / m(cgs)^2 * 4*pi
        cnst *= accel;
        double dn = dn_b/1.0e6; // dn - CGS cm^-3



        double ee_ev = en_a*J_2_EV; // test ptl's kinetic energy in [eV]

        // calculate psi(x), f, g, gp
        double dumb = vprt/vt; // MKS/MKS  ,  v/v_th
        double dd = dumb*dumb; // dd = x, dumb = sqrt(x)
        double d3 = dumb*dd;
        double dum = d3*(ap0 + ap1*dd + ap2*dd*dd + ap3*dd*dd*dd + ap4*dd*dd*dd*dd);
        double ap_psi = 1.0 - 1.0/(1.0 + dum); //psi(x) = Phi(v/vth)

        double f = (2.0 - 1.0/dd)*ap_psi+2.257*dumb*exp(-dd);
        double g = 2.0*mass*ap_psi/mass_b;
        double gp = 2.257*exp(-dd)*(mass/mass_b)*dumb;

        // find Coulomb logarithm
        //  bmax  = 7.4d-3*dsqrt(en_b_ev/1000d0*1.d13/(charge**2*dn_b))
        //  bmin1 = 1.4d-10*charge*charge_b/(ee_ev + en_b_ev)*1000d0
        //  massab = mass_au*mass_b/(mass_au+mass_b)
        //  bmin2 = 6.2d-10/(dsqrt(massab*1836.d0*en_b_ev/1000d0))


        double bmax = 7.4e-3*sqrt(en_b_ev*1.e10/(charge*charge*dn));
        double bmin1 = 1.4e-7*charge*charge_b/(ee_ev+en_b_ev);
        double massab = mass*mass_b/(mass + mass_b);
        double bmin2 = 6.2e-10/(sqrt(massab*1836.*en_b_ev/1.e3));
        double bmin = max(bmin1,bmin2);
        double clog = log(bmax/bmin);

        // collision frequencies - scattering, slowing down
        double dnu_b = cnst*dn*charge_b*charge_b;
        freq_scat = abs(clog*dnu_b*f);
        freq_slow = abs(clog*dnu_b*g);
        freq_fac0 = en_a*(1.0 - mass_b*gp/(g*mass));
    }

    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{
        if(use_varying_bg && vb.is_in_range(magnetic_field, r, z, psi)){
            vb.interp_bg_profile(psi,theta,species.nonadiabatic_idx, den, temp, up);
        } else {
            den = species.eq_den.value(magnetic_field, psi, r, z);
            temp = species.eq_temp.value(magnetic_field, psi, r, z);
            up = species.eq_flow.value(magnetic_field, psi, r, z);

            int ftype=species.eq_flow_type;
            up=(ftype>=1)? up*r : up;
            if(ftype>=2) { // ftype==2 might slow down the performance a litte.
                double phi=0; // assign random angle assuming toroidal symmetry.
                double br,bz,bphi;
                magnetic_field.bvec_interpol(r,z,phi,br,bz,bphi);
                up *= bphi/sqrt(br*br+bz*bz+bphi*bphi);
            }
        }
    }

    public:

    MonteCarloCollider(){}

    MonteCarloCollider(NLReader::NamelistReader& nlr, const MagneticField<DeviceType>& magnetic_field, int n_nonadiabatic_species){
        nlr.use_namelist("col_param");
        period = nlr.get<int>("col_period", 3);
        accel = nlr.get<bool>("col_accel", false);

        accel_n = nlr.get<int>("col_accel_n",2);
        accel_factor1 = nlr.get<double>("col_accel_factor1",10.0);
        accel_factor2 = nlr.get<double>("col_accel_factor2",10.0);

        accel_pin1 = nlr.get<double>("col_accel_pin1", magnetic_field.inpsi);
        double psi_range = magnetic_field.outpsi - magnetic_field.inpsi;
        accel_pout1 = nlr.get<double>("col_accel_pout1", magnetic_field.inpsi + 0.1*psi_range);
        accel_pin2 = nlr.get<double>("col_accel_pin2", magnetic_field.outpsi - 0.1*psi_range);
        accel_pout2 = nlr.get<double>("col_accel_pout2", magnetic_field.outpsi);

        // Normalize
        accel_pin1  *= magnetic_field.equil.xpt_psi;
        accel_pout1 *= magnetic_field.equil.xpt_psi;
        accel_pin2  *= magnetic_field.equil.xpt_psi;
        accel_pout2 *= magnetic_field.equil.xpt_psi;

        en_col_on = (nlr.get<int>("col_en_col_on", 1))==1;
        moving_frame = nlr.get<bool>("col_moving_frame", true);

        use_varying_bg = (nlr.get<int>("col_varying_bg", 0))==1;
        if(use_varying_bg){
            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");

            // Initialize varying background
            vb = VaryingBackground<Device>(nlr, magnetic_field, n_nonadiabatic_species);

            if(vb.period<period) exit_XGC("\nError: must set col_vb_period >= col_period\n");
        }
    }

    void update_vb_all_species(const MagneticField<DeviceType>& magnetic_field, Plasma& plasma) {
        plasma.for_all_nonadiabatic_species([&](Species<DeviceType>& species){
            // Copy particles to the device
            TIMER("copy_ptl_to_device",
                species.copy_particles_to_device_if_not_resident() );

            // A snapshot is created for ions and electrons every col_vb_period calls to collision
            TIMER("COL_SNAPSHOT",
                vb.update(magnetic_field, species) );
        });
    }
};

void monte_carlo_collisions(int istep, const MagneticField<DeviceType>& magnetic_field, Plasma& plasma, MonteCarloCollider<DeviceType>& col_mc, double dt_in);

#endif