XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
analytic_psi.hpp
Go to the documentation of this file.
1 #ifndef ANALYTIC_PSI_HPP
2 #define ANALYTIC_PSI_HPP
3 
4 #include "globals.hpp"
5 #include "equil.hpp"
6 
7 // Provides setup and get psi for an analytic example of psi
8 // get_analytic_psi() will be replaced with bicubic interp eventually
9 
10 namespace{
11 
12 double get_circular_analytic_psi(const Equilibrium& equil, double r, double z){
13  double rp = (r - equil.axis.r)/equil.bounds.get_r_width();
14  double zp = (z - equil.axis.z)/equil.bounds.get_z_width();
15 
16  return rp*rp + zp*zp;
17 }
18 
19 
20 double get_analytic_psi(const Equilibrium& equil, double r, double z){
21  double rp = r - equil.axis.r;
22  double zp = z - equil.axis.z;
23  double r_xpt1p = equil.xpt.r-equil.axis.r;
24  double z_xpt1p = equil.xpt.z-equil.axis.z;
25  double r_xpt2p = equil.xpt2.r-equil.axis.r;
26  double z_xpt2p = equil.xpt2.z-equil.axis.z;
27  double zx1p = z - equil.xpt.z;
28  double zx2p = z - equil.xpt2.z;
29  double avg_z_xpt = (equil.xpt2.z - equil.xpt.z)/2;
30  double avg_psi_xpt = (equil.xpt2_psi + equil.xpt_psi)/2;
31  double ar = 1/pow(avg_z_xpt,2);
32 
33  double r_xptp = (zp>0 ? r_xpt2p : r_xpt1p);
34  double z_xptp = (zp>0 ? z_xpt2p : z_xpt1p);
35  double zxp = (zp>0 ? zx2p : zx1p);
36  double xpt_psi = (zp>0 ? equil.xpt2_psi : equil.xpt_psi);
37 
38  double a = -2/pow(z_xptp,3);
39  double b = -3/pow(z_xptp,2);
40  if (zp>z_xpt2p || zp<z_xpt1p){
41  a = 0;
42  b = -1/pow(z_xptp,2);
43  }
44 
45 
46  // Setup if no second x-point
47  double no_xpt_psi = 0;
48  double anx = 0;
49  if ((zp>0 && !equil.set_xpt2)){
50  // Use xpoint 1 as reference for scale
51  xpt_psi = equil.xpt_psi;
52  no_xpt_psi = equil.xpt_psi;
53  anx = 1/pow(z_xpt1p,2);
54  // Remove xpt coefficients:
55  a = 0;
56  b = 0;
57  }
58 
59  return ar*pow(rp-r_xptp/z_xptp*zp,2)*avg_psi_xpt + (a*pow(zxp,3) + b*pow(zxp,2))*xpt_psi + xpt_psi
60  + anx*pow(zp,2)*no_xpt_psi - no_xpt_psi;
61 }
62 
63 void get_analytic_psi_derivs(const Equilibrium& equil, double r, double z, double& dpsidr, double& dpsidz){
64  double rp = r - equil.axis.r;
65  double zp = z - equil.axis.z;
66  double r_xpt1p = equil.xpt.r-equil.axis.r;
67  double z_xpt1p = equil.xpt.z-equil.axis.z;
68  double r_xpt2p = equil.xpt2.r-equil.axis.r;
69  double z_xpt2p = equil.xpt2.z-equil.axis.z;
70  double zx1p = z - equil.xpt.z;
71  double zx2p = z - equil.xpt2.z;
72  double avg_z_xpt = (equil.xpt2.z - equil.xpt.z)/2;
73  double avg_psi_xpt = (equil.xpt2_psi + equil.xpt_psi)/2;
74  double ar = 1/pow(avg_z_xpt,2);
75 
76  double r_xptp = (zp>0 ? r_xpt2p : r_xpt1p);
77  double z_xptp = (zp>0 ? z_xpt2p : z_xpt1p);
78  double zxp = (zp>0 ? zx2p : zx1p);
79  double xpt_psi = (zp>0 ? equil.xpt2_psi : equil.xpt_psi);
80 
81  double a = -2/pow(z_xptp,3);
82  double b = -3/pow(z_xptp,2);
83  if (zp>z_xpt2p || zp<z_xpt1p){
84  a = 0;
85  b = -1/pow(z_xptp,2);
86  }
87 
88 
89  // Setup if no second x-point
90  double no_xpt_psi = 0;
91  double anx = 0;
92  if ((zp>0 && !equil.set_xpt2)){
93  // Use xpoint 1 as reference for scale
94  xpt_psi = equil.xpt_psi;
95  no_xpt_psi = equil.xpt_psi;
96  anx = 1/pow(z_xpt1p,2);
97  // Remove xpt coefficients:
98  a = 0;
99  b = 0;
100  }
101 
102  dpsidr = 2*ar*(rp-r_xptp/z_xptp*zp)*avg_psi_xpt;
103  dpsidz = -2*ar*(rp-r_xptp/z_xptp*zp)*avg_psi_xpt*r_xptp/z_xptp + (3*a*pow(zxp,2) + 2*b*zxp)*xpt_psi + 2*anx*zp*no_xpt_psi;
104 }
105 
106 
107 void get_analytic_psi_second_derivs(const Equilibrium& equil, double z, double& dpsidrdz, double& dpsidr2, double& dpsidz2){
108  double zp = z - equil.axis.z;
109  double r_xpt1p = equil.xpt.r-equil.axis.r;
110  double z_xpt1p = equil.xpt.z-equil.axis.z;
111  double r_xpt2p = equil.xpt2.r-equil.axis.r;
112  double z_xpt2p = equil.xpt2.z-equil.axis.z;
113  double zx1p = z - equil.xpt.z;
114  double zx2p = z - equil.xpt2.z;
115  double avg_z_xpt = (equil.xpt2.z - equil.xpt.z)/2;
116  double avg_psi_xpt = (equil.xpt2_psi + equil.xpt_psi)/2;
117  double ar = 1/pow(avg_z_xpt,2);
118 
119  double r_xptp = (zp>0 ? r_xpt2p : r_xpt1p);
120  double z_xptp = (zp>0 ? z_xpt2p : z_xpt1p);
121  double zxp = (zp>0 ? zx2p : zx1p);
122  double xpt_psi = (zp>0 ? equil.xpt2_psi : equil.xpt_psi);
123 
124  double a = -2/pow(z_xptp,3);
125  double b = -3/pow(z_xptp,2);
126  if (zp>z_xpt2p || zp<z_xpt1p){
127  a = 0;
128  b = -1/pow(z_xptp,2);
129  }
130 
131 
132  // Setup if no second x-point
133  double no_xpt_psi = 0;
134  double anx = 0;
135  if ((zp>0 && !equil.set_xpt2)){
136  // Use xpoint 1 as reference for scale
137  xpt_psi = equil.xpt_psi;
138  no_xpt_psi = equil.xpt_psi;
139  anx = 1/pow(z_xpt1p,2);
140  // Remove xpt coefficients:
141  a = 0;
142  b = 0;
143  }
144 
145  dpsidrdz = -2*ar*(r_xptp/z_xptp)*avg_psi_xpt;
146  dpsidr2 = 2*ar*avg_psi_xpt;
147  dpsidz2 = 2*ar*(r_xptp/z_xptp)*avg_psi_xpt*r_xptp/z_xptp + (6*a*zxp + 2*b)*xpt_psi + 2*anx*no_xpt_psi;
148 }
149 
150 void plot_equilibrium(Equilibrium& equil){
151  int nx = 80;
152  int ny = 40;
153  double dr = equil.bounds.get_r_width()/nx;
154  double dz = equil.bounds.get_z_width()/ny;
155  printf("\n");
156  std::vector<int> prevrow(nx+1,-1);
157 
158  for (int j=ny;j>=0;j--){
159  printf("\n");
160  int prevnum = -1;
161  for (int i=0;i<=nx;i++){
162  double r = equil.bounds.min_r + i*dr;
163  double z = equil.bounds.min_z + j*dz;
164  int num = get_analytic_psi(equil,r,z)*5.0;
165  char mychar = 48+num;
166  if (prevnum!=num || prevrow[i]!=num) printf("%c",mychar);
167  else printf(" ");
168  prevnum = num;
169  prevrow[i] = num;
170  }
171  }
172  printf("\n");
173 }
174 
175 // Create a psi eq grid based off the analytic psi
176 Kokkos::View<double**, Kokkos::LayoutRight, HostType> analytic_psi_grid(const Equilibrium& equil,
177  const View<double*, HostType>& r,
178  const View<double*, HostType>& z){
179 
180  Kokkos::View<double**, Kokkos::LayoutRight, HostType> psi_grid("psi_grid",z.size(),r.size());
181  for (int j=0; j<z.size(); j++){
182  for (int i=0; i<r.size(); i++){
183  psi_grid(j,i) = get_analytic_psi(equil, r(i), z(j));
184  }
185  }
186  return psi_grid;
187 }
188 
189 // Another psi eq grid for testing: Simple, circular case
190 Kokkos::View<double**, Kokkos::LayoutRight, HostType> circular_analytic_psi_grid(const Equilibrium& equil,
191  const View<double*, HostType>& r,
192  const View<double*, HostType>& z){
193 
194  Kokkos::View<double**, Kokkos::LayoutRight, HostType> psi_grid("psi_grid",z.size(),r.size());
195  for (int j=0; j<z.size(); j++){
196  for (int i=0; i<r.size(); i++){
197  psi_grid(j,i) = get_circular_analytic_psi(equil, r(i), z(j));
198  }
199  }
200  return psi_grid;
201 }
202 
203 }
204 #endif
Equilibrium::xpt2
RZPair xpt2
coordinates of 2nd X-point
Definition: equil.hpp:87
RZPair::z
double z
Definition: grid_structs.hpp:30
Equilibrium::set_xpt2
bool set_xpt2
Whether to use a 2nd X-point.
Definition: equil.hpp:86
RZBounds::min_r
double min_r
Definition: rz_bounds.hpp:5
Equilibrium::bounds
RZBounds bounds
Min and max for r and z.
Definition: equil.hpp:80
Equilibrium::xpt2_psi
double xpt2_psi
psi coordinate at 2nd X-point
Definition: equil.hpp:89
Equilibrium::xpt
RZPair xpt
coordinates of 1st X-point
Definition: equil.hpp:84
RZPair::r
double r
Definition: grid_structs.hpp:29
RZBounds::get_z_width
KOKKOS_INLINE_FUNCTION double get_z_width() const
Definition: rz_bounds.hpp:28
RZBounds::get_r_width
KOKKOS_INLINE_FUNCTION double get_r_width() const
Definition: rz_bounds.hpp:24
RZBounds::min_z
double min_z
Definition: rz_bounds.hpp:7
Equilibrium::axis
RZPair axis
Definition: equil.hpp:90
Equilibrium
Definition: equil.hpp:29
Equilibrium::xpt_psi
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:81
Generated by   doxygen 1.8.5