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 mu, double c2_2m){
#ifdef XGC_FLAT_B
    double B=eq_axis_b;
#else
    double br, bz, bphi;
    magnetic_field.bvec_interpol(r,z,0.0,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 SimdGyroRadius;

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

    Simd<double> rho[2];

    public:

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

    // calculate rho in constructor
    template<class Device>
    KOKKOS_INLINE_FUNCTION SimdGyroRadius(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector2D &x, Simd<double>& phi, Simd<double>& mu, Simd<int>& itr){
        // Need to follow fields to planes rather than midpoint
        // follow_field takes a vector for phi_dest (since electrons could have different ones)
        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[i_simd])*grid.delta_phi;
            phimax[i_simd] = phimin[i_simd] + grid.delta_phi;
        }
        SimdVector2D xp[2];
        magnetic_field.follow_field(x,phi,phimin,xp[0]);
        magnetic_field.follow_field(x,phi,phimax,xp[1]);
        Simd<double> psip[2];
        magnetic_field.get_psi(xp[0],psip[0]);
        magnetic_field.get_psi(xp[1],psip[1]);

        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            for (int iphi=0; iphi<=1; iphi++){
                double r = xp[iphi].r[i_simd];
                double z = xp[iphi].z[i_simd];
                double ti=species.eq_temp.value(magnetic_field, psip[iphi][i_simd],r,z);
                double br, bz, bphi;
                magnetic_field.bvec_interpol(r,z,0.0,br,bz,bphi);
                double B = sqrt( br*br + bz*bz + bphi*bphi );
                rho[iphi][i_simd]=sqrt(2.0*B*mu[i_simd]/(ti*EV_2_J));
            }
        }
    }

};

#else
template<>
struct SimdGyroRadius<GyroKin>{
    private:

    Simd<double> rho;

    public:

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

    // calculate rho in constructor
    template<class Device>
    KOKKOS_INLINE_FUNCTION SimdGyroRadius(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector2D &x, Simd<double>& phi, Simd<double>& mu, Simd<int>& itr){
        for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
            rho[i_simd]=gyro_radius(magnetic_field,x.r[i_simd],x.z[i_simd],mu[i_simd],species.c2_2m);
        }
    }
};
#endif

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

    // Default constructor:
    KOKKOS_INLINE_FUNCTION SimdGyroRadius(){}

    // Empty constructor to look like the GyroKin constructor:
    template<class Device>
    KOKKOS_INLINE_FUNCTION SimdGyroRadius(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector2D &x, Simd<double>& phi, Simd<double>& mu, Simd<int>& itr){ }
};


#endif