XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
gyro_radius.hpp
Go to the documentation of this file.
1 #ifndef GYRO_RADIUS_HPP
2 #define GYRO_RADIUS_HPP
3 
4 #include "globals.hpp"
5 #include "basic_physics.hpp"
6 #include "magnetic_field.hpp"
7 #include "grid.hpp"
8 #include "species.hpp"
9 
10 template<class Device>
11 KOKKOS_INLINE_FUNCTION double gyro_radius(const MagneticField<Device> &magnetic_field, double r, double z, double mu, double c2_2m){
12 #ifdef XGC_FLAT_B
13  double B=eq_axis_b;
14 #else
15  double br, bz, bphi;
16  magnetic_field.bvec_interpol(r,z,0.0,br,bz,bphi);
17  double B = sqrt( br*br + bz*bz + bphi*bphi );
18 #endif
19  return gyro_radius(B, mu, c2_2m);
20 }
21 
22 // Generic declaration
23 template<KinType PT>
25 
26 #if defined(XGC1) && defined(NEWGYROMATRIX)
27 template<>
28 struct SimdGyroRadius<GyroKin>{
29  private:
30 
31  Simd<double> rho[2];
32 
33  public:
34 
35  // Allows array-like access pattern
36  KOKKOS_INLINE_FUNCTION Simd<double> operator [](int i) const {return rho[i];}
37  KOKKOS_INLINE_FUNCTION Simd<double>& operator [](int i) {return rho[i];}
38 
39  // calculate rho in constructor
40  template<class Device>
41  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){
42  // Need to follow fields to planes rather than midpoint
43  // follow_field takes a vector for phi_dest (since electrons could have different ones)
44  Simd<double> phimin;
45  Simd<double> phimax;
46  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
47  phimin[i_simd] = grid.get_plane_index(phi[i_simd])*grid.delta_phi;
48  phimax[i_simd] = phimin[i_simd] + grid.delta_phi;
49  }
50  SimdVector2D xp[2];
51  magnetic_field.follow_field(x,phi,phimin,xp[0]);
52  magnetic_field.follow_field(x,phi,phimax,xp[1]);
53  Simd<double> psip[2];
54  magnetic_field.get_psi(xp[0],psip[0]);
55  magnetic_field.get_psi(xp[1],psip[1]);
56 
57  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
58  for (int iphi=0; iphi<=1; iphi++){
59  double r = xp[iphi].r[i_simd];
60  double z = xp[iphi].z[i_simd];
61  double ti=species.eq_temp.value(magnetic_field, psip[iphi][i_simd],r,z);
62  double br, bz, bphi;
63  magnetic_field.bvec_interpol(r,z,0.0,br,bz,bphi);
64  double B = sqrt( br*br + bz*bz + bphi*bphi );
65  rho[iphi][i_simd]=sqrt(2.0*B*mu[i_simd]/(ti*EV_2_J));
66  }
67  }
68  }
69 
70 };
71 
72 #else
73 template<>
75  private:
76 
78 
79  public:
80 
81  // Allows array-like access pattern
82  KOKKOS_INLINE_FUNCTION double operator [](int i) const {return rho[i];}
83  KOKKOS_INLINE_FUNCTION double& operator [](int i) {return rho[i];}
84 
85  // calculate rho in constructor
86  template<class Device>
87  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){
88  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
89  rho[i_simd]=gyro_radius(magnetic_field,x.r[i_simd],x.z[i_simd],mu[i_simd],species.c2_2m);
90  }
91  }
92 };
93 #endif
94 
95 template<>
97  // Empty for Drift Kinetic (no gyroradius)
98 
99  // Default constructor:
100  KOKKOS_INLINE_FUNCTION SimdGyroRadius(){}
101 
102  // Empty constructor to look like the GyroKin constructor:
103  template<class Device>
104  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){ }
105 };
106 
107 
108 #endif
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:229
KOKKOS_INLINE_FUNCTION void follow_field(const SimdVector2D &x_org, const Simd< double > &phi_org, const Simd< double > &phi_dest, SimdVector2D &x_dest) const
Definition: magnetic_field.tpp:295
Definition: gyro_radius.hpp:24
double c2_2m
c2/2m
Definition: species.hpp:87
constexpr double EV_2_J
Conversion rate ev to J.
Definition: constants.hpp:5
Definition: globals.hpp:89
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)
Definition: gyro_radius.hpp:87
Definition: magnetic_field.hpp:12
real(kind=8) function gyro_radius(x, mu, isp)
Definition: poisson_extra.F90:74
Definition: grid.hpp:67
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector2D &x, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:82
KOKKOS_INLINE_FUNCTION SimdGyroRadius()
Definition: gyro_radius.hpp:100
Simd< double > z
Definition: simd.hpp:141
Definition: globals.hpp:90
KOKKOS_INLINE_FUNCTION int get_plane_index(double phi) const
Definition: grid.tpp:872
Simd< double > r
Definition: simd.hpp:140
Definition: magnetic_field.F90:1
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)
Definition: gyro_radius.hpp:104
Definition: simd.hpp:139
Definition: species.hpp:75
Simd< double > rho
Definition: gyro_radius.hpp:77
Eq::Profile< Device > eq_temp
Definition: species.hpp:124
double delta_phi
Distance between planes.
Definition: grid.hpp:271