XGCa
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
profile.hpp
Go to the documentation of this file.
1 #ifndef PROFILE_HPP
2 #define PROFILE_HPP
3 
4 #include "magnetic_field.hpp"
5 
6 // Putting this in a namespace since Shape names sound so generic
7 namespace Eq {
8 
9 enum Shape{
19 };
20 
21 //parameters for linear equal distance interpolation of profiles
22 //maybe need some compiler flags for generality. Currently hard coded to 1000 and true.
23 const int ftn_lin_n=1000; // 1000 data point for linear interpolation
24 
25 // data structure for linear interpolation -- used in eq_ftn_type
26 template<class Device>
28  int n;
29  double p_min,p_del;
30  Kokkos::View<double*,Kokkos::LayoutRight,Device> v;
31 
32  public:
33 
34  CustomLinShape(double p_min_in, double p_max_in, double* v_in)
35  : n(ftn_lin_n),
36  p_min(p_min_in),
37  p_del((p_max_in - p_min_in)/(ftn_lin_n-1)),
38  v("v",ftn_lin_n)
39  {
40  array_deep_copy(v,v_in);
41  }
42 
44 
45  KOKKOS_INLINE_FUNCTION double interp(double x) const{
46  double xn=(x-p_min)/p_del;
47  xn=max(min(xn,double(n)),0.0);
48  int ip=int(xn);
49  ip=min(ip,n-2);
50  double wp=1.0-(xn-double(ip));
51  return v(ip)*wp+v(ip+1)*(1.0-wp);
52  }
53 
54  KOKKOS_INLINE_FUNCTION double dinterp(double x) const{
55  double xn=(x-p_min)/p_del;
56  xn=max(min(xn,double(n)),0.0);
57  int ip=int(xn);
58  ip=min(ip,n-2);
59  return (v(ip+1)-v(ip))/p_del;
60  }
61 };
62 
63 // data structure for equilbirum profile
64 template<class Device>
65 class Profile{
67  double inx[5];
68  double iny[4];
69  double sv[6];
70  double p_min;
71  double p_max;
72 
73  // for linear interpolation
75 
76  // get value of profile at the given coordinates
77  KOKKOS_INLINE_FUNCTION double value_at_coords(const MagneticField<DeviceType> &b_field, double psi,double r,double z) const {
78  if(shape==Constant){
79  return iny[0];
80  } else if(shape==TanH){
81  return sv[1]*tanh((inx[0]-psi)*sv[2]) + sv[0];
82  } else if(shape==Linear){
83  return sv[0]*psi+sv[1];
84  } else if(shape==Shape5){
85  if(abs(psi/b_field.equil.xpt_psi)<1.0e-3) psi = 1.0e-3*b_field.equil.xpt_psi;
86  double tmp2=sv[2]*(inx[0]-sqrt(inx[0]*psi)); // z
87  double tmp3=exp(tmp2); // expz
88  double tmp4=exp(-tmp2); // expmz
89  // A * ( (1+z*slope)*expz - expmz )/(expz + expmz) + B
90  return sv[1]*((1+tmp2*sv[3])*tmp3-tmp4)/(tmp3+tmp4) + sv[0];
91  } else if(shape==CustomLinear){
92  psi=min(max(psi,p_min),p_max);
93  return lin.interp(psi);
94  } else{ // Shouldnt get here
95  return 0.0;
96  }
97  }
98 
99  public:
100 
101  // Profile from Fortran
102  Profile(double* inx_in, double* iny_in, double* sv_in, double p_min_in, double p_max_in, int shape_int, double* v_in )
103  : p_min(p_min_in),
104  p_max(p_max_in)
105  {
106  for (int i=0; i<5; i++) inx[i]=inx_in[i];
107  for (int i=0; i<4; i++) iny[i]=iny_in[i];
108  for (int i=0; i<6; i++) sv[i]=sv_in[i];
109 
110  // Currently, only a few shapes available
111  if(shape_int==0){
112  shape=Constant;
113  } else if(shape_int==1){
114  shape=TanH;
115  } else if(shape_int==2){
116  shape=Linear;
117  } else if(shape_int==5){
118  shape=Shape5;
119  } else if(shape_int==-1){ // shape==-1 is also CustomSpline; this is caught in the init right now
122  } else {
123  printf("invalid shape number in eq_ftn_setup: %d",shape_int);
124  exit(0);
125  }
126  }
127 
128  // Constructor for simple linear profile
129  Profile(double val_at_zero_psi, double slope)
130  : shape(Linear)
131  {
132  for (int i=0; i<5; i++) inx[i]=0.0;
133  for (int i=0; i<4; i++) iny[i]=0.0;
134  for (int i=0; i<6; i++) sv[i]=0.0;
135  sv[0] = slope;
136  sv[1] = val_at_zero_psi;
137  }
138 
139  // Profile for analytic testing
140  Profile(Shape shape_in)
141  : shape(shape_in),
142  p_min(0.0),
143  p_max(0.0)
144  {
145  for (int i=0; i<5; i++) inx[i]=0.0;
146  for (int i=0; i<4; i++) iny[i]=0.0;
147  for (int i=0; i<6; i++) sv[i]=0.0;
148  if(shape==Constant){
149  iny[0] = 2.0; // value
150  } else if(shape==TanH){
151  inx[0] = 1.0;
152  sv[0] = 3.0; // half height
153  sv[1] = 3.0; // value at psi=inx[0]
154  sv[2] = 5.0; // width
155  } else if(shape==Linear){
156  sv[0] = 2.0; // slope
157  sv[1] = 3.0; // value at psi=0
158  } else if(shape==Shape5){
159  inx[0] = 2.5e-1;
160  sv[0] = 5.5e2;
161  sv[1] = 4.5e2;
162  sv[2] = 2.1;
163  sv[3] = 8.0e-3;
164  } else if(shape==CustomLinear){
165  double v_in[ftn_lin_n];
166  p_min = 0.0;
167  p_max = 1.0;
168  for (int i = 0; i<ftn_lin_n; i++){
169  v_in[i] = 1.0+i/(ftn_lin_n-1.0);
170  }
172  }
173  }
174 
176 
177  // get value of profile, including decay if out of range
178  KOKKOS_INLINE_FUNCTION double value(const MagneticField<DeviceType> &b_field, double psi_in,double r,double z) const {
179  bool is_in_rgn12 = b_field.equil.is_in_region_1_or_2(r,z,psi_in);
180  if(is_in_rgn12 && (psi_in<=b_field.outpsi)){
181  return value_at_coords(b_field,psi_in,r,z);
182  } else if(is_in_rgn12 && (psi_in>b_field.outpsi)){
183  // Get function value at psi=b_field.outpsi
184  // let ftn decay exponentially to eq_out_decay_factor
185  // times the value at b_field.outpsi with a decay
186  // length of eq_out_decay_width
187  double psi=b_field.outpsi;
188  double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
189  double tmp3=tmp2*(1.0-b_field.equil.out_decay_factor);
190  return tmp3*exp(-abs(psi_in-psi)/b_field.equil.out_decay_width)+b_field.equil.out_decay_factor*tmp2;
191  } else if(!is_in_rgn12 && z<b_field.equil.xpt_z){
192  // For primary private flux region
193  double psi=b_field.equil.xpt_psi;
194  double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
195  double tmp3=tmp2*(1.0-b_field.equil.priv_flux_decay_factor);
196  return tmp3*exp(-abs(psi_in-psi)/b_field.equil.priv_flux_decay_width)+b_field.equil.priv_flux_decay_factor*tmp2;
197  } else if(!is_in_rgn12 && z>b_field.equil.xpt2_z && b_field.equil.set_xpt2){
198  // For primary private flux region
199  double psi=b_field.equil.xpt2_psi;
200  double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
201  double tmp3=tmp2*(1.0-b_field.equil.priv_flux_decay_factor);
202  return tmp3*exp(-abs(psi_in-psi)/b_field.equil.priv_flux_decay_width)+b_field.equil.priv_flux_decay_factor*tmp2;
203  } else {
204  // No way to get here? - ALS
205  double psi=b_field.equil.xpt_psi;
206  return value_at_coords(b_field,psi,r,z);
207  }
208  }
209 
210  // get derivative of profile
211  KOKKOS_INLINE_FUNCTION double slope(const MagneticField<DeviceType> &b_field, double psi_in,double r,double z) const {
212  // region searching and enforcing region 3 value to x-point value
213  //rh return derivative=0 for psi values larger than sml_outpsi
214  //rh This is not a serious constraint since eq_dftn is called only in
215  //rh the particle push to get the turbulence drive for delta-f simulations
216  if(!b_field.equil.is_in_region_1_or_2(r,z,psi_in) || (psi_in>b_field.outpsi)){
217  return 0.0;
218  }
219 
220  if(shape==Constant){
221  return 0.0;
222  } else if(shape==TanH){
223  double tmp2 = 1.0/cosh((inx[0]-psi_in)*sv[2]);
224  return -sv[2]*sv[1]*tmp2*tmp2;
225  } else if(shape==Linear){
226  return sv[0];
227  } else if(shape==Shape5){
228  double tmp2=sv[2]*(inx[0]-sqrt(inx[0]*psi_in)); // z
229  double tmp3=exp(tmp2); // expz
230  // A * ( (1+z*slope)*expz - expmz )/(expz + expmz) + B
231  if(abs(psi_in/b_field.equil.xpt_psi)<1.0e-3){
232  return 0.0;
233  } else {
234  return -sv[1]*tmp2*tmp3*tmp3*(4.0+sv[3]*(1.0+tmp3*tmp3+2.0*tmp2))
235  /(2.0*(1.0+tmp2*tmp2)*(1.0+tmp2*tmp2)*(sqrt(inx[0]*psi_in)-psi_in));
236  }
237  } else if(shape==CustomLinear){
238  if(p_min<psi_in && psi_in<p_max){
239  return lin.dinterp(psi_in);
240  } else {
241  return 0.0;
242  }
243  } else{ // Shouldnt get here
244  return 0.0;
245  }
246  }
247 
248 };
249 }
250 #endif
const int ftn_lin_n
Definition: profile.hpp:23
double p_del
Definition: profile.hpp:29
KOKKOS_INLINE_FUNCTION double slope(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:211
KOKKOS_INLINE_FUNCTION double interp(double x) const
Definition: profile.hpp:45
Shape shape
Definition: profile.hpp:66
Definition: profile.hpp:11
Definition: magnetic_field.hpp:9
Profile(double *inx_in, double *iny_in, double *sv_in, double p_min_in, double p_max_in, int shape_int, double *v_in)
Definition: profile.hpp:102
Definition: profile.hpp:12
double p_min
Definition: profile.hpp:70
int n
Definition: profile.hpp:28
CustomLinShape< Device > lin
Definition: profile.hpp:74
CustomLinShape(double p_min_in, double p_max_in, double *v_in)
Definition: profile.hpp:34
Definition: profile.hpp:15
double inx[5]
Definition: profile.hpp:67
Definition: profile.hpp:14
double outpsi
Boundary condition used in a few spots.
Definition: magnetic_field.hpp:26
KOKKOS_INLINE_FUNCTION double value_at_coords(const MagneticField< DeviceType > &b_field, double psi, double r, double z) const
Definition: profile.hpp:77
Definition: profile.hpp:17
KOKKOS_INLINE_FUNCTION double dinterp(double x) const
Definition: profile.hpp:54
Definition: profile.hpp:13
Definition: profile.hpp:16
Profile()
Definition: profile.hpp:175
double p_min
Definition: profile.hpp:29
double iny[4]
Definition: profile.hpp:68
CustomLinShape()
Definition: profile.hpp:43
Equilibrium< Device > equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
Shape
Definition: profile.hpp:9
Profile(double val_at_zero_psi, double slope)
Definition: profile.hpp:129
Definition: profile.hpp:18
KOKKOS_INLINE_FUNCTION double value(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:178
Definition: profile.hpp:27
Profile(Shape shape_in)
Definition: profile.hpp:140
Definition: profile.hpp:10
Kokkos::View< double *, Kokkos::LayoutRight, Device > v
Definition: profile.hpp:30
double sv[6]
Definition: profile.hpp:69
double p_max
Definition: profile.hpp:71
Definition: profile.hpp:65
void array_deep_copy(T *array, const Kokkos::View< T *, Kokkos::LayoutRight, DeviceType > &view)
Definition: array_deep_copy.hpp:11