col_species.hpp

Go to the documentation of this file.

#ifndef COL_SPECIES_HPP
#define COL_SPECIES_HPP

#include <vector>
#include <Kokkos_Core.hpp>
#include <Kokkos_DualView.hpp>

#include "vgrid_distribution.hpp"

extern "C" double get_den_at_node(int node, int species);
extern "C" double get_t_ev_at_node(int node, int species);

struct CollisionSpeciesScalars{
    int gid;               // Corresponding XGC species
    int grid_index;        // Which grid this species is on
    bool is_electron;      // Whether species is electron
    double T_avg;
    double mass;
    double mass_au;
    double charge_eu;
    double den_lambda_gamma;
    double conv_factor;
    double numeric_vth2;            //[col_fm_core_init<col_fm_core]
    double numeric_t;               //[col_fm_core_init<col_fm_core]
    double dens;
    double ens;
    double mom;                     // [col_fm_core_init<col_fm_core]

    //> Computes the base collision frequency -- Eq. (6) in Hager et al.
    KOKKOS_INLINE_FUNCTION double lambda_gamma_pair(const CollisionSpeciesScalars& sp_b ) const {// C(a->b)
        double lambda;
        // NRL Plasma Formulary 2011 version, p.34
        //(c) Mixed ion-ion collisions (here, same species ion-ion collisions)
        if(is_electron || sp_b.is_electron){
            //if electron is invovled
            if(is_electron && sp_b.is_electron){
                //(a) Thermal electron-electron collisions
                // Note that pointers for colliding and target are same.
                lambda = 2.35e1 - log( sqrt(den_lambda_gamma*1e-6)*pow(T_avg,-1.25))
                                - sqrt(std::abs(1e-5 + pow(log(T_avg)-2,2)/16));
            } else {
                double ti_ev,massi,densi,chargei,te_ev,masse,dense;
                if(is_electron){
                    ti_ev = sp_b.T_avg;
                    massi = sp_b.mass_au;
                    densi = sp_b.den_lambda_gamma;
                    chargei = sp_b.charge_eu;
                    te_ev = T_avg;
                    masse = mass_au;
                    dense = den_lambda_gamma;
                } else {
                    te_ev = sp_b.T_avg;
                    masse = sp_b.mass_au;
                    dense = sp_b.den_lambda_gamma;
                    ti_ev = T_avg;
                    massi = mass_au;
                    densi = den_lambda_gamma;
                    chargei = charge_eu;
                }
                //(b) Electron-ion (and ion-electron) collisions
                if(ti_ev*masse/massi<te_ev){
                    if(te_ev<1e1*(chargei*chargei)){
                        lambda = 2.3e1 - log(chargei*sqrt(dense*1e-6/(te_ev*te_ev*te_ev)));
                    } else {
                        lambda = 2.4e1 - log(sqrt(dense*1e-6)/te_ev);
                    }
                } else {
                    lambda = 3e1-log((chargei*chargei)*sqrt(densi*1e-6/(ti_ev*ti_ev*ti_ev))/massi);
                }
            }
        } else {
            lambda = 2.3e1 - log(charge_eu*sp_b.charge_eu*(mass_au+sp_b.mass_au)
                                      /(mass_au*sp_b.T_avg+sp_b.mass_au*T_avg)
                                      *sqrt( den_lambda_gamma*1e-6*(charge_eu*charge_eu)/T_avg
                                            +sp_b.den_lambda_gamma*1e-6*(sp_b.charge_eu*sp_b.charge_eu)/sp_b.T_avg));
        }
        double tmp = lambda*(charge_eu*charge_eu)*(sp_b.charge_eu*sp_b.charge_eu)*pow(UNIT_CHARGE,4)/(pow(8.8542e-12,2)*8*PI);
        return tmp/mass; // colliding : a, target : b
    }
};

template<class Device>
class CollisionSpecies{
    public:
    const VGridDistribution<HostType> f;

    Kokkos::DualView<CollisionSpeciesScalars**,Device> s;

    Kokkos::DualView<double****,Kokkos::LayoutRight,Device> pdf_n;    //[get_local_f<setup_col_sp] Given initial PDF
    Kokkos::DualView<double****,Kokkos::LayoutRight,Device> pdf_np1;  //[get_local_f<setup_col_sp] To be used as PDF at next Picard step. BUT has a role of "df" at end

    CollisionSpecies<Device>(){}

    CollisionSpecies<Device>(const VGridDistribution<HostType>& f_in, std::vector<Species<Device>> &all_species, int mb_n_nodes)
        : f(f_in),
          s("s", f.n_species(), mb_n_nodes),
          pdf_n("pdf_n", mb_n_nodes, f.n_species(), f.n_vz(), f.n_vr()),
          pdf_np1("pdf_np1", mb_n_nodes, f.n_species(), f.n_vz(), f.n_vr())
    {
        int col_species_cnt=0;
        for (int isp = 0; isp<all_species.size(); isp++){
            if (!all_species[isp].is_adiabatic){
                for (int mesh_ind = 0; mesh_ind<mb_n_nodes; mesh_ind++){
                    s.h_view(col_species_cnt, mesh_ind).gid = isp;
                    s.h_view(col_species_cnt, mesh_ind).is_electron = all_species[isp].is_electron;
                    s.h_view(col_species_cnt, mesh_ind).grid_index = all_species[isp].collision_grid_index;
                }
                col_species_cnt++;
            }
        }
    }

    //> Convert the distribution function to the units used
    // internally by the collision operator
    // The distribution function set up with f0_module is in normalized velocity space
    // whereas the collision operator uses SI units.
    void setup_one(int isp, int mesh_ind, const Species<Device> &species, int local_node_ind, int first_node) {
        s.h_view(isp,mesh_ind).mass = species.mass;
        s.h_view(isp,mesh_ind).mass_au = s.h_view(isp,mesh_ind).mass/PROTON_MASS;
        s.h_view(isp,mesh_ind).charge_eu = species.charge/UNIT_CHARGE;

        double conv_inv_mass=TWOPI*UNIT_CHARGE/s.h_view(isp,mesh_ind).mass;
        s.h_view(isp,mesh_ind).den_lambda_gamma = get_den_at_node(local_node_ind + first_node, s.h_view(isp,mesh_ind).gid);
        s.h_view(isp,mesh_ind).T_avg = get_t_ev_at_node(local_node_ind + first_node, s.h_view(isp,mesh_ind).gid);
        s.h_view(isp,mesh_ind).conv_factor = 1.0/sqrt(s.h_view(isp,mesh_ind).T_avg*(conv_inv_mass*conv_inv_mass*conv_inv_mass));

        // Local f with basic correction:
        // Simply cut off all negative values
        for (int imu = 0; imu < pdf_n.h_view.extent(3); imu++){
            double smu_n = f.get_smu_n(imu);
            for (int invp = 0; invp<pdf_n.h_view.extent(2); invp++){
                pdf_n.h_view(mesh_ind,isp,invp,imu) = std::max(f(species.nonadiabatic_idx, imu,local_node_ind,invp),0.0)*s.h_view(isp,mesh_ind).conv_factor/smu_n;
            }
        }

        // pdf_np1  =  pdf_n - replace w deep copy - ALS
        for (int iz = 0; iz < pdf_np1.h_view.extent(2); iz++){
            for (int ir = 0; ir < pdf_np1.h_view.extent(3); ir++){
                pdf_np1.h_view(mesh_ind,isp,iz,ir) = pdf_n.h_view(mesh_ind,isp,iz,ir);
            }
        }
    }

    void setup_all(const std::vector<Species<Device>>& all_species, Kokkos::View<int**,HostType>& mesh_nodes, int mesh_batch_ind, int first_node){

        for (int mesh_ind = 0; mesh_ind<s.extent(1); mesh_ind++){
            for (int isp = 0; isp < s.extent(0); isp++){
                setup_one(isp, mesh_ind, all_species[s.h_view(isp,mesh_ind).gid], mesh_nodes(mesh_batch_ind,mesh_ind), first_node);
            }
        }

#ifdef USE_GPU
        // Copy dual views to device
        Kokkos::deep_copy(s.d_view, s.h_view);
        Kokkos::deep_copy(pdf_n.d_view, pdf_n.h_view);
#endif
    }

    void lambda_gamma_pair(double dt, Kokkos::View<double***,Device>& gammac){
        // Done on host, so create a mirror
        //typename Kokkos::View<double***,Device>::HostMirror 
        auto gammac_h = Kokkos::create_mirror_view(gammac);

        for (int mesh_ind = 0; mesh_ind<s.extent(1); mesh_ind++){
            // gamma
            for (int isp = 0; isp < s.extent(0); isp++){
                for (int jsp = 0; jsp < s.extent(0); jsp++){
                    gammac_h(mesh_ind,jsp,isp) = dt*s.h_view(isp,mesh_ind).lambda_gamma_pair(s.h_view(jsp,mesh_ind));
                }
            }
        }
        Kokkos::deep_copy(gammac, gammac_h);
    }

    inline int n() const{
        return f.n_species();
    }
};

#endif