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 SimdGyroWeights<GyroKin>{
29  //private:
30 
31  SimdLinearWeights rho_wts[2];
32 
33  //public:
34 
35  // Allows array-like access pattern
36  //KOKKOS_INLINE_FUNCTION Simd<double> operator [](int i) const {return rho_wts[i];}
37  //KOKKOS_INLINE_FUNCTION Simd<double>& operator [](int i) {return rho_wts[i];}
38 
39  // calculate rho in constructor
40  template<class Device>
41  KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector &v, Simd<double>& mu){
42 
43  // Need to follow fields to planes rather than midpoint
44  // follow_field takes a vector for phi_dest (since electrons could have different ones)
45  // Get toroidal angle
46  Simd<double> phi_mod = v.phi;
47  grid.wedge_modulo_phi(phi_mod);
48  Simd<double> phimin;
49  Simd<double> phimax;
50  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
51  phimin[i_simd] = grid.get_plane_index(phi_mod[i_simd])*grid.delta_phi;
52  phimax[i_simd] = phimin[i_simd] + grid.delta_phi;
53  }
54  SimdVector2D xp[2];
55  magnetic_field.follow_field(v.x(),phi_mod,phimin,xp[0]);
56  magnetic_field.follow_field(v.x(),phi_mod,phimax,xp[1]);
57 
58  for (int iphi=0; iphi<2; iphi++){
59  Simd<double> psip;
60  magnetic_field.get_psi(xp[iphi], iphi=0 ? phimin : phimax, psip);
61  Simd<double> rho;
62  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
63  double r = xp[iphi].r[i_simd];
64  double z = xp[iphi].z[i_simd];
65  double ti=species.eq_temp.value(magnetic_field, psip[i_simd],r,z);
66  double br, bz, bphi;
67  magnetic_field.bvec_interpol(r,z,0.0,br,bz,bphi);
68  double B = sqrt( br*br + bz*bz + bphi*bphi );
69  rho[i_simd]=sqrt(2.0*B*mu[i_simd]/(ti*EV_2_J));
70  }
71 
72  rho_wts[i].set(rho, species.gyro_avg_matrices.inv_drho, species.gyro_avg_matrices.nrho);
73  }
74  }
75 
76  KOKKOS_INLINE_FUNCTION int irho(int i_simd) const{
77  // If NEWGYROMATRIX, there are 2 rho wts so must specify (even though the i indices are identical)
78  return rho_wts[0].i[i_simd];
79  }
80 };
81 
82 #else
83 template<>
85  //private:
86 
88 
89  //public:
90 
91  // Allows array-like access pattern
92  //KOKKOS_INLINE_FUNCTION double operator [](int i) const {return rho_wts[i];}
93  //KOKKOS_INLINE_FUNCTION double& operator [](int i) {return rho_wts[i];}
94 
95  // calculate rho in constructor
96  template<class Device>
97  KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector &v, Simd<double>& mu){
98  Simd<double> rho;
99 
100  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
101  rho[i_simd]=gyro_radius(magnetic_field,v.r[i_simd],v.z[i_simd],mu[i_simd],species.c2_2m);
102  }
103 
104  rho_wts.set(rho, species.gyro_avg_matrices.inv_drho, species.gyro_avg_matrices.nrho);
105  }
106 
107  KOKKOS_INLINE_FUNCTION int irho(int i_simd) const{
108  return rho_wts.i[i_simd];
109  }
110 };
111 #endif
112 
113 template<>
115  // Empty for Drift Kinetic (no gyroradius)
116 
117  // Default constructor:
118  KOKKOS_INLINE_FUNCTION SimdGyroWeights(){}
119 
120  // Empty constructor to look like the GyroKin constructor:
121  template<class Device>
122  KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid<Device> &grid, const MagneticField<Device> &magnetic_field, const Species<DeviceType> &species, const SimdVector &v, Simd<double>& mu){ }
123 
124  KOKKOS_INLINE_FUNCTION int irho(int i_simd) const{
125  return 0;
126  }
127 };
128 
129 
130 #endif
Simd< double > r
Definition: simd.hpp:150
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:321
Definition: simd.hpp:149
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector &v, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:105
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:403
Definition: linear_weights.hpp:46
double c2_2m
c2/2m
Definition: species.hpp:87
KOKKOS_INLINE_FUNCTION int irho(int i_simd) const
Definition: gyro_radius.hpp:107
constexpr double EV_2_J
Conversion rate ev to J.
Definition: constants.hpp:5
Definition: globals.hpp:89
Definition: magnetic_field.hpp:13
real(kind=8) function gyro_radius(x, mu, isp)
Definition: poisson_extra.F90:74
KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid< Device > &grid, const MagneticField< Device > &magnetic_field, const Species< DeviceType > &species, const SimdVector &v, Simd< double > &mu)
Definition: gyro_radius.hpp:97
Definition: grid.hpp:21
KOKKOS_INLINE_FUNCTION SimdGyroWeights(const Grid< Device > &grid, const MagneticField< Device > &magnetic_field, const Species< DeviceType > &species, const SimdVector &v, Simd< double > &mu)
Definition: gyro_radius.hpp:122
int nrho
Definition: gyro_avg_mat.hpp:16
KOKKOS_INLINE_FUNCTION SimdVector2D & x()
Definition: simd.hpp:156
SimdLinearWeights rho_wts
Definition: gyro_radius.hpp:87
KOKKOS_INLINE_FUNCTION SimdGyroWeights()
Definition: gyro_radius.hpp:118
Definition: gyro_radius.hpp:24
Simd< double > z
Definition: simd.hpp:141
Definition: globals.hpp:90
double inv_drho
Definition: gyro_avg_mat.hpp:18
KOKKOS_INLINE_FUNCTION int get_plane_index(double phi) const
Definition: grid.tpp:133
KOKKOS_INLINE_FUNCTION void wedge_modulo_phi(Simd< double > &phi_mod) const
Definition: grid.tpp:86
Simd< double > phi
Definition: simd.hpp:152
Simd< double > z
Definition: simd.hpp:151
Simd< double > r
Definition: simd.hpp:140
Definition: magnetic_field.F90:1
Definition: simd.hpp:139
KOKKOS_INLINE_FUNCTION double value(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:319
GyroAverageMatrices< HostType > gyro_avg_matrices
Definition: species.hpp:147
Definition: species.hpp:75
Eq::Profile< Device > eq_temp
Definition: species.hpp:131
double delta_phi
Distance between planes.
Definition: grid.hpp:157
KOKKOS_INLINE_FUNCTION int irho(int i_simd) const
Definition: gyro_radius.hpp:124