gyro_radius.hpp

Go to the documentation of this file.

#ifndef GYRO_RADIUS_HPP
#define GYRO_RADIUS_HPP

#include "globals.hpp"
#include "basic_physics.hpp"
#include "magnetic_field.hpp"
#include "grid.hpp"
#include "species.hpp"

template<class Device>
KOKKOS_INLINE_FUNCTION double gyro_radius(const MagneticField<Device> &magnetic_field, double r, double z, double phi, double mu, double c2_2m){
#ifdef XGC_FLAT_B
    double B=eq_axis_b;
#else
    double br, bz, bphi;
    magnetic_field.bvec_interpol(r,z,phi,br,bz,bphi);
    double B = sqrt( br*br + bz*bz + bphi*bphi );
#endif
    return gyro_radius(B, mu, c2_2m);
}

// Generic declaration
template<KinType PT>
struct SimdGyroWeights;

#if defined(XGC1) && defined(NEWGYROMATRIX)
template<>
struct SimdGyroWeights<GyroKin>{
    //private:

    SimdLinearWeights rho_wts[2];

    //public:

    // Allows array-like access pattern
    //KOKKOS_INLINE_FUNCTION Simd<double> operator [](int i) const {return rho_wts[i];}
    //KOKKOS_INLINE_FUNCTION Simd<double>& operator [](int i) {return rho_wts[i];}

    // calculate rho in constructor
    template<class Device>
    KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector &v, Simd<double>& mu){

        // Need to follow fields to planes rather than midpoint
        // follow_field takes a vector for phi_dest (since electrons could have different ones)
        // Get toroidal angle
        Simd<double> phi_mod = v.phi;
        grid.wedge_modulo_phi(phi_mod);
        Simd<double> phimin;
        Simd<double> phimax;
        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            phimin[i_simd] = grid.get_plane_index(phi_mod[i_simd])*grid.delta_phi;
            phimax[i_simd] = phimin[i_simd] + grid.delta_phi;
        }
        SimdVector2D xp[2];
        magnetic_field.follow_field(v.x(),phi_mod,phimin,xp[0]);
        magnetic_field.follow_field(v.x(),phi_mod,phimax,xp[1]);

        for (int iphi=0; iphi<2; iphi++){
            Simd<double> psip;
            magnetic_field.get_psi(xp[iphi], iphi=0 ? phimin : phimax, psip);
            Simd<double> rho;
            for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
                double r = xp[iphi].r[i_simd];
                double z = xp[iphi].z[i_simd];
                double ti=species.eq_temp.value(magnetic_field, psip[i_simd],r,z);
                double br, bz, bphi;
                magnetic_field.bvec_interpol(r,z,phi_mod[i_simd],br,bz,bphi);
                double B = sqrt( br*br + bz*bz + bphi*bphi );
                rho[i_simd]=sqrt(2.0*B*mu[i_simd]/(ti*EV_2_J));
            }

            rho_wts[i].set(rho, species.gyro_avg_matrices.inv_drho, species.gyro_avg_matrices.nrho);
        }
    }

    KOKKOS_INLINE_FUNCTION int irho(int i_simd) const{
        // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
        return rho_wts[0].i[i_simd];
    }
};

#else
template<>
struct SimdGyroWeights<GyroKin>{
    //private:

    SimdLinearWeights rho_wts;

    //public:

    // Allows array-like access pattern
    //KOKKOS_INLINE_FUNCTION double operator [](int i) const {return rho_wts[i];}
    //KOKKOS_INLINE_FUNCTION double& operator [](int i) {return rho_wts[i];}

    KOKKOS_INLINE_FUNCTION SimdGyroWeights(){}

    // calculate rho in constructor
    template<class Device>
    KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector &v, Simd<double>& mu){
        Simd<double> rho;

        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            rho[i_simd]=gyro_radius(magnetic_field,v.r[i_simd],v.z[i_simd],v.phi[i_simd],mu[i_simd],species.c2_2m);
        }

        rho_wts.set(rho, species.gyro_avg_matrices.inv_drho, species.gyro_avg_matrices.nrho);
    }

    KOKKOS_INLINE_FUNCTION int irho(int i_simd) const{
        return rho_wts.i[i_simd];
    }
};
#endif

template<>
struct SimdGyroWeights<DriftKin>{
    // Empty for Drift Kinetic (no gyroradius)

    // Default constructor:
    KOKKOS_INLINE_FUNCTION SimdGyroWeights(){}

    // Empty constructor to look like the GyroKin constructor:
    template<class Device>
    KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector &v, Simd<double>& mu){ }

    KOKKOS_INLINE_FUNCTION int irho(int i_simd) const{
        return 0;
    }
};


#endif