XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator 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 "magnetic_field.hpp"
6 #include "grid.hpp"
7 #include "species.hpp"
8 
9 template<class Device>
10 KOKKOS_INLINE_FUNCTION double gyro_radius(const MagneticField<Device> &magnetic_field, double r, double z, double mu, double c2_2m){
11 #ifdef XGC_FLAT_B
12  double B=eq_axis_b;
13 #else
14  double br, bz, bphi;
15  magnetic_field.bvec_interpol(r,z,0.0,br,bz,bphi);
16  double B = sqrt( br*br + bz*bz + bphi*bphi );
17 #endif
18  return sqrt(mu/B/c2_2m);
19 }
20 
21 // Generic declaration
22 template<KinType PT>
24 
25 #if defined(XGC1) && defined(NEWGYROMATRIX)
26 template<>
27 struct SimdGyroRadius<GyroKin>{
28  private:
29 
30  Simd<double> rho[2];
31 
32  public:
33 
34  // Allows array-like access pattern
35  KOKKOS_INLINE_FUNCTION double operator [](int i) const {return rho[i];}
36  KOKKOS_INLINE_FUNCTION double& operator [](int i) {return rho[i];}
37 
38  // calculate rho in constructor
39  template<class Device>
40  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){
41  // Need to follow fields to planes rather than midpoint
42  // follow_field takes a vector for phi_dest (since electrons could have different ones)
43  Simd<double> phimin;
44  Simd<double> phimax;
45  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
46  phimin[i_simd] = grid.phimin;
47  phimax[i_simd] = grid.phimax;
48  }
49  SimdVector2D xp[2];
50  magnetic_field.follow_field(x,phi,phimin,xp[0]);
51  magnetic_field.follow_field(x,phi,phimax,xp[1]);
52  Simd<double> psip[2];
53  magnetic_field.get_psi(xp[0].r,xp[0].z,psip[0]);
54  magnetic_field.get_psi(xp[1].r,xp[1].z,psip[1]);
55 
56  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
57  for (int iphi=0; iphi<=1; iphi++){
58  double r = xp[iphi].r[i_simd];
59  double z = xp[iphi].z[i_simd];
60  double ti=species.eq_temp.value(magnetic_field, psip[iphi][i_simd],r,z);
61  double br, bz, bphi;
62  magnetic_field.bvec_interpol(r,z,0.0,br,bz,bphi);
63  double B = sqrt( br*br + bz*bz + bphi*bphi );
64  rho[iphi][i_simd]=sqrt(2.0*B*mu[i_simd]/(ti*EV_2_J));
65  }
66  }
67  }
68 
69 };
70 
71 #else
72 template<>
74  private:
75 
77 
78  public:
79 
80  // Allows array-like access pattern
81  KOKKOS_INLINE_FUNCTION double operator [](int i) const {return rho[i];}
82  KOKKOS_INLINE_FUNCTION double& operator [](int i) {return rho[i];}
83 
84  // calculate rho in constructor
85  template<class Device>
86  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){
87  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
88  rho[i_simd]=gyro_radius(magnetic_field,x.r[i_simd],x.z[i_simd],mu[i_simd],species.c2_2m);
89  }
90  }
91 };
92 #endif
93 
94 template<>
96  // Empty for Drift Kinetic (no gyroradius)
97 
98  // Default constructor:
99  KOKKOS_INLINE_FUNCTION SimdGyroRadius(){}
100 
101  // Empty constructor to look like the GyroKin constructor:
102  template<class Device>
103  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){ }
104 };
105 
106 
107 #ifdef NEWGYROMATRIX
108 // With NEWGYROMATRIX, save the results.
109 template<class Device>
110 KOKKOS_INLINE_FUNCTION void save_rhon(const Grid<Device> &grid,SimdGyroRadius<GyroKin>& rho,int i_item,const TmpSpecies<Device> &tmp_species){
111  for (int i_simd = 0; i_simd<SIMD_SIZE; i_simd++){
112  // Index for saving rhon for ion push
113  int i = i_item*SIMD_SIZE + i_simd;
114 # ifdef XGC1
115  tmp_species.rhon(i,0) = double(grid.nrho)*min(rho[0][i_simd]/grid.rhomax,1.0);
116  tmp_species.rhon(i,1) = double(grid.nrho)*min(rho[1][i_simd]/grid.rhomax,1.0);
117 # else
118  tmp_species.rhon(i,0) = double(grid.nrho)*min(rho[i_simd]/grid.rhomax,1.0);
119 # endif
120  }
121 }
122 #endif
123 
124 #endif
KOKKOS_INLINE_FUNCTION void bvec_interpol(double r, double z, double phi, double &br, double &bz, double &bphi) const
Definition: magnetic_field.tpp:150
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:204
Definition: gyro_radius.hpp:23
double c2_2m
c2/2m
Definition: species.hpp:24
KOKKOS_INLINE_FUNCTION void get_psi(const Simd< double > &r_in, const Simd< double > &z_in, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:37
const double EV_2_J
Conversion rate ev to J.
Definition: globals.hpp:36
Definition: globals.hpp:14
Definition: magnetic_field.hpp:9
Definition: grid.hpp:10
int nrho
Definition: grid.hpp:127
double phimin
Definition: grid.hpp:130
KOKKOS_INLINE_FUNCTION SimdGyroRadius()
Definition: gyro_radius.hpp:99
real(kind=8) function gyro_radius(x, mu, isp)
Definition: poisson_extra_common.F90:74
Simd< double > z
Definition: globals.hpp:56
Definition: globals.hpp:15
double phimax
Definition: grid.hpp:131
double rhomax
Definition: grid.hpp:128
Simd< double > r
Definition: globals.hpp:55
Definition: magnetic_field.F90:1
Definition: globals.hpp:54
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)
Definition: gyro_radius.hpp:86
Definition: species.hpp:13
Definition: gyro_radius.hpp:73
Simd< double > rho
Definition: gyro_radius.hpp:76
Eq::Profile< Device > eq_temp
Definition: species.hpp:39
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)
Definition: gyro_radius.hpp:103
Definition: species.hpp:71