XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
profile.hpp
Go to the documentation of this file.
1 #ifndef PROFILE_HPP
2 #define PROFILE_HPP
3 
4 #include "file_reader.hpp"
5 #include "magnetic_field.hpp"
6 #ifdef USE_MPI
7 #include "my_mpi.hpp"
8 #endif
9 
10 extern "C" void init_ftn_lin(int num, double* psi, double* var, double psi_min, double dpsi, int n, double* v);
11 
12 // Putting this in a namespace since Shape names sound so generic
13 namespace Eq {
14 
15 enum Shape{
25 };
26 
27 // data structure for linear interpolation
28 template<class Device>
30  static constexpr int n = 1000; // 1000 data point for linear interpolation
31  double p_min;
32  double p_max;
33  double p_del;
34  Kokkos::View<double*,Kokkos::LayoutRight,Device> v;
35 
36  public:
37 
38  CustomLinShape<Device>(double psi_normalization, const std::string& filename)
39  : v(NoInit("v"),n)
40  {
41  // Only one MPI rank reads the files, then broadcasts to the others
42  int num;
43  View<double*, CLayout, HostType> psi;
44  View<double*, CLayout, HostType> var;
45  if(is_rank_zero()){
46  // Open file
47  FileReader file_reader;
48  file_reader.open(filename);
49  file_reader.read(&num, 1);
50  if(num<=0) exit_XGC("\nError in profile ftn init : invalid number of data\n");
51 
52  psi = View<double*, CLayout, HostType>(NoInit("psi"), num);
53  var = View<double*, CLayout, HostType>(NoInit("var"), num);
54 
55  // Read data
56  for(int i=0; i<num; i++){
57  file_reader.read(&psi(i), 1);
58  file_reader.read(&var(i), 1);
59  }
60 
61  // Check end of file
62  int flag;
63  file_reader.read(&flag, 1);
64  if(flag!=-1) exit_XGC("\nError in profile ftn init : invalid number of data or ending flag -1 is not set correctly\n");
65  }
66 
67  // Broadcast file size to other processes
68 #ifdef USE_MPI
69  int root_rank = 0;
70  MPI_Bcast(&num, 1, MPI_INT, root_rank, SML_COMM_WORLD);
71 #endif
72  if(!is_rank_zero()){
73  psi = View<double*, CLayout, HostType>(NoInit("psi"), num);
74  var = View<double*, CLayout, HostType>(NoInit("var"), num);
75  }
76  // Broadcast arrays
77 #ifdef USE_MPI
78  MPI_Bcast( psi.data(), psi.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
79  MPI_Bcast( var.data(), var.size(), MPI_DOUBLE, root_rank, SML_COMM_WORLD);
80 #endif
81 
82  // Normalization
83  for(int i=0; i<num; i++){
84  psi(i) *= psi_normalization;
85  }
86 
87  // save psi range
88  p_min = psi(0);
89  p_max = psi(num-1);
90  p_del = (p_max - p_min)/(n-1);
91 
92  // Set value of v on host
93  View<double*,CLayout,HostType> v_h(NoInit("v_h"), v.layout());
94 #ifdef PSPLINE
95  init_ftn_lin(num, psi.data(), var.data(), p_min, p_del, n, v_h.data());
96 #else
97  exit_XGC("\nError: CustomLinear profile shape requires PSPLINE.\n");
98 #endif
99 
100  // Copy to device
101  Kokkos::deep_copy(v, v_h);
102  }
103 
104  // Constructor provided an array of values and psis
105  CustomLinShape<Device>(const View<double*, CLayout, HostType>& psi, const View<double*, CLayout, HostType>& var)
106  : v(NoInit("v"),n)
107  {
108  int num = var.size();
109 
110  // save psi range
111  p_min = psi(0);
112  p_max = psi(num-1);
113  p_del = (p_max - p_min)/(n-1);
114 
115  // Set value of v on host
116  View<double*,CLayout,HostType> v_h(NoInit("v_h"), v.layout());
117 #ifdef PSPLINE
118  init_ftn_lin(num, psi.data(), var.data(), p_min, p_del, n, v_h.data());
119 #else
120  exit_XGC("\nError: CustomLinear profile shape requires PSPLINE.\n");
121 #endif
122 
123  // Copy to device
124  Kokkos::deep_copy(v, v_h);
125  }
126 
127  // Constructor with linear value for testing
128  CustomLinShape<Device>(double p_min_in, double p_max_in)
129  : v(NoInit("v"),n),
130  p_min(p_min_in),
131  p_max(p_max_in),
132  p_del((p_max_in-p_min_in)/(n-1))
133  {
134  View<double*,CLayout,HostType> v_h(NoInit("v_h"), v.layout());
135  for (int i = 0; i<n; i++){
136  v_h(i) = 1.0 + i/(n-1.0);
137  }
138 
139  // Copy to device
140  Kokkos::deep_copy(v, v_h);
141  }
142 
144 
145  KOKKOS_INLINE_FUNCTION double interp(double x) const{
146  // Bound by psi limits - this seems redundant with the bounds immediately following
147  x = min(max(x,p_min),p_max);
148  double xn=(x-p_min)/p_del;
149  xn=max(min(xn,double(n)),0.0);
150  int ip=int(xn);
151  ip=min(ip,n-2);
152  double wp=1.0-(xn-double(ip));
153  return v(ip)*wp+v(ip+1)*(1.0-wp);
154  }
155 
156  KOKKOS_INLINE_FUNCTION double dinterp(double x) const{
157  if(p_min<x && x<p_max){
158  double xn=(x-p_min)/p_del;
159  xn=max(min(xn,double(n)),0.0);
160  int ip=int(xn);
161  ip=min(ip,n-2);
162  return (v(ip+1)-v(ip))/p_del;
163  } else {
164  return 0.0;
165  }
166  }
167 };
168 
169 // Data structure for equilibrium profile
170 template<class Device>
171 class Profile{
172  public:
173  static constexpr int N_X = 5;
174  static constexpr int N_Y = 4;
175  static constexpr int N_V = 6;
176 
177  private:
178 
180  double inx[N_X];
181  double iny[N_Y];
182  double sv[N_V];
183 
184  // for linear interpolation
186 
187  // get value of profile at the given coordinates
188  KOKKOS_INLINE_FUNCTION double value_at_coords(const MagneticField<DeviceType> &b_field, double psi,double r,double z) const {
189  if(shape==Constant){
190  return iny[0];
191  } else if(shape==TanH){
192  return sv[1]*tanh((inx[0]-psi)*sv[2]) + sv[0];
193  } else if(shape==Linear){
194  return sv[0]*psi+sv[1];
195  } else if(shape==Shape5){
196  if(abs(psi/b_field.equil.xpt_psi)<1.0e-3) psi = 1.0e-3*b_field.equil.xpt_psi;
197  double tmp2=sv[2]*(inx[0]-sqrt(inx[0]*psi)); // z
198  double tmp3=exp(tmp2); // expz
199  double tmp4=exp(-tmp2); // expmz
200  // A * ( (1+z*slope)*expz - expmz )/(expz + expmz) + B
201  return sv[1]*((1+tmp2*sv[3])*tmp3-tmp4)/(tmp3+tmp4) + sv[0];
202  } else if(shape==CustomLinear){
203  return lin.interp(psi);
204  } else{ // Shouldnt get here
205  return 0.0;
206  }
207  }
208 
209  public:
210 
211  // This function converts from the integer currently provided in the namelist
212  // into the enum
213  static Shape shape_from_int(int shape_int){
214  if(shape_int==0){
215  return Constant;
216  } else if(shape_int==1){
217  return TanH;
218  } else if(shape_int==2){
219  return Linear;
220  } else if(shape_int==5){
221  return Shape5;
222  } else if(shape_int==-1){ // shape==-1 is also CustomSpline; this is caught in the init right now
223  return CustomLinear;
224  } else {
225  printf("Invalid shape number in profile setup: %d",shape_int);
226  exit_XGC("\nError: Invalid shape number\n");
227  return Constant; // Not reached
228  }
229  }
230 
231  Profile(double psi_normalization, Shape shape_in, const std::vector<double>& inx_in, const std::vector<double>& iny_in, const std::string& filename)
232  : shape(shape_in)
233  {
234  // Copy vector inputs into member C-style arrays
235  for(int i=0; i<inx_in.size(); i++) inx[i] = inx_in[i];
236  for(int i=0; i<iny_in.size(); i++) iny[i] = iny_in[i];
237 
238 
239  // Precalculations depending on profile shape
240  if(shape==Constant){
241  // No additional setup needed
242  }else if(shape==TanH){
243  sv[0]=0.5*(iny[0]+iny[1]);
244  sv[1]=0.5*(iny[0]-iny[1]);
245  sv[2]=2/inx[1];
246  }else if(shape==Linear){
247  sv[0]=(iny[1]-iny[0])/ inx[1]; // slope
248  sv[1]=iny[0]-sv[0]*inx[0]; // offset
249  }else if(shape==Shape5){
250  // A ( (1+z slope) e^z - e^-z )/(e^z + e^-z) + B
251  // z= 4 ( center - sqrt(center psi) )/ width
252  // iny(1) - edge_val : ped_top
253  // iny(2) - out_val : ped bottom
254  // inx(1) - center , inx(2) - width
255  // iny(3) - core(axis)_val , inx(3) - core(axis)_psi
256  sv[0]=0.5*(iny[0]+iny[1]);
257  sv[1]=0.5*(iny[0]-iny[1]); // height
258  sv[2]=4/inx[1]; // width inverse
259  sv[3]=(iny[2]-iny[0]) / (sqrt(inx[2]) - sqrt(inx[0]))
260  /( - (sv[1]+1.0e-99) * sv[2] * sqrt(inx[0]) ); // slope
261  }else if(shape==CustomLinear){
262  lin = CustomLinShape<Device>(psi_normalization, filename);
263  }else{
264  exit_XGC("\nError: Invalid profile shape provided.\n");
265  }
266  }
267 
268  // Constructor that takes a pre-existing custom profile x and var
269  Profile(const View<double*, CLayout, HostType>& psi, const View<double*, CLayout, HostType>& var)
270  : shape(CustomLinear),
271  lin(CustomLinShape<Device>(psi, var))
272  {}
273 
274  // Constructor for simple linear profile
275  Profile(double val_at_zero_psi, double slope)
276  : shape(Linear)
277  {
278  for (int i=0; i<5; i++) inx[i]=0.0;
279  for (int i=0; i<4; i++) iny[i]=0.0;
280  for (int i=0; i<6; i++) sv[i]=0.0;
281  sv[0] = slope;
282  sv[1] = val_at_zero_psi;
283  }
284 
285  // Profile for analytic testing
286  Profile(Shape shape_in)
287  : shape(shape_in)
288  {
289  for (int i=0; i<5; i++) inx[i]=0.0;
290  for (int i=0; i<4; i++) iny[i]=0.0;
291  for (int i=0; i<6; i++) sv[i]=0.0;
292  if(shape==Constant){
293  iny[0] = 2.0; // value
294  } else if(shape==TanH){
295  inx[0] = 1.0;
296  sv[0] = 3.0; // half height
297  sv[1] = 3.0; // value at psi=inx[0]
298  sv[2] = 5.0; // width
299  } else if(shape==Linear){
300  sv[0] = 2.0; // slope
301  sv[1] = 3.0; // value at psi=0
302  } else if(shape==Shape5){
303  inx[0] = 2.5e-1;
304  sv[0] = 5.5e2;
305  sv[1] = 4.5e2;
306  sv[2] = 2.1;
307  sv[3] = 8.0e-3;
308  } else if(shape==CustomLinear){
309  // Linear value for testing
310  double p_min = 0.0;
311  double p_max = 1.0;
312  lin = CustomLinShape<Device>(p_min, p_max);
313  }
314  }
315 
317 
318  // get value of profile, including decay if out of range
319  KOKKOS_INLINE_FUNCTION double value(const MagneticField<DeviceType> &b_field, double psi_in,double r,double z) const {
320  bool is_in_rgn12 = b_field.equil.is_in_region_1_or_2(r,z,psi_in);
321  if(is_in_rgn12 && (psi_in<=b_field.outpsi)){
322  return value_at_coords(b_field,psi_in,r,z);
323  } else if(is_in_rgn12 && (psi_in>b_field.outpsi)){
324  // Get function value at psi=b_field.outpsi
325  // let ftn decay exponentially to eq_out_decay_factor
326  // times the value at b_field.outpsi with a decay
327  // length of eq_out_decay_width
328  double psi=b_field.outpsi;
329  double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
330  double tmp3=tmp2*(1.0-b_field.equil.out_decay_factor);
331  return tmp3*exp(-abs(psi_in-psi)/b_field.equil.out_decay_width)+b_field.equil.out_decay_factor*tmp2;
332  } else if(!is_in_rgn12 && z<b_field.equil.xpt_z){
333  // For primary private flux region
334  double psi=b_field.equil.xpt_psi;
335  double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
336  double tmp3=tmp2*(1.0-b_field.equil.priv_flux_decay_factor);
337  return tmp3*exp(-abs(psi_in-psi)/b_field.equil.priv_flux_decay_width)+b_field.equil.priv_flux_decay_factor*tmp2;
338  } else if(!is_in_rgn12 && z>b_field.equil.xpt2_z && b_field.equil.set_xpt2){
339  // For primary private flux region
340  double psi=b_field.equil.xpt2_psi;
341  double tmp2=value_at_coords(b_field,psi,b_field.equil.axis_r,b_field.equil.axis_z);
342  double tmp3=tmp2*(1.0-b_field.equil.priv_flux_decay_factor);
343  return tmp3*exp(-abs(psi_in-psi)/b_field.equil.priv_flux_decay_width)+b_field.equil.priv_flux_decay_factor*tmp2;
344  } else {
345  // No way to get here? - ALS
346  double psi=b_field.equil.xpt_psi;
347  return value_at_coords(b_field,psi,r,z);
348  }
349  }
350 
351  // get derivative of profile
352  KOKKOS_INLINE_FUNCTION double slope(const MagneticField<DeviceType> &b_field, double psi_in,double r,double z) const {
353  // region searching and enforcing region 3 value to x-point value
354  //rh return derivative=0 for psi values larger than sml_outpsi
355  //rh This is not a serious constraint since eq_dftn is called only in
356  //rh the particle push to get the turbulence drive for delta-f simulations
357  if(!b_field.equil.is_in_region_1_or_2(r,z,psi_in) || (psi_in>b_field.outpsi)){
358  return 0.0;
359  }
360 
361  if(shape==Constant){
362  return 0.0;
363  } else if(shape==TanH){
364  double tmp2 = 1.0/cosh((inx[0]-psi_in)*sv[2]);
365  return -sv[2]*sv[1]*tmp2*tmp2;
366  } else if(shape==Linear){
367  return sv[0];
368  } else if(shape==Shape5){
369  double tmp2=sv[2]*(inx[0]-sqrt(inx[0]*psi_in)); // z
370  double tmp3=exp(tmp2); // expz
371  // A * ( (1+z*slope)*expz - expmz )/(expz + expmz) + B
372  if(abs(psi_in/b_field.equil.xpt_psi)<1.0e-3){
373  return 0.0;
374  } else {
375  return -inx[0]*sv[1]*sv[2]*tmp3*tmp3*(4.0+sv[3]*(1.0+tmp3*tmp3+2.0*tmp2))
376  /(2.0*(1.0+tmp3*tmp3)*(1.0+tmp3*tmp3)*(sqrt(inx[0]*psi_in)));
377  }
378  } else if(shape==CustomLinear){
379  return lin.dinterp(psi_in);
380  } else{ // Shouldnt get here
381  return 0.0;
382  }
383  }
384 
385 };
386 }
387 #endif
double xpt2_z
z coordinate at 2nd X-point
Definition: equil.hpp:87
double p_del
Definition: profile.hpp:33
static constexpr int N_X
Definition: profile.hpp:173
bool is_rank_zero()
Definition: globals.hpp:27
MPI_Comm SML_COMM_WORLD
Definition: my_mpi.cpp:4
KOKKOS_INLINE_FUNCTION double slope(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:352
KOKKOS_INLINE_FUNCTION double interp(double x) const
Definition: profile.hpp:145
double iny[N_Y]
Definition: profile.hpp:181
Shape shape
Definition: profile.hpp:179
Definition: profile.hpp:17
static constexpr int N_V
Definition: profile.hpp:175
Definition: magnetic_field.hpp:12
double out_decay_width
width for exponential decay for psi&gt;sml_outpsi
Definition: equil.hpp:95
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
Definition: profile.hpp:18
void read(T *var, int n)
Definition: file_reader.hpp:28
static Shape shape_from_int(int shape_int)
Definition: profile.hpp:213
CustomLinShape< Device > lin
Definition: profile.hpp:185
double axis_z
z coordinate of axis
Definition: equil.hpp:90
bool set_xpt2
Whether to use a 2nd X-point.
Definition: equil.hpp:84
Definition: profile.hpp:21
Definition: profile.hpp:20
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:188
double sv[N_V]
Definition: profile.hpp:182
double out_decay_factor
profiles decay exponentially to f(sml_outpsi)/decay_factor for psi&gt;sml_outpsi
Definition: equil.hpp:93
Definition: profile.hpp:23
KOKKOS_INLINE_FUNCTION double dinterp(double x) const
Definition: profile.hpp:156
Definition: profile.hpp:19
double axis_r
r coordinate of axis
Definition: equil.hpp:89
Profile(double psi_normalization, Shape shape_in, const std::vector< double > &inx_in, const std::vector< double > &iny_in, const std::string &filename)
Definition: profile.hpp:231
Definition: profile.hpp:22
Profile()
Definition: profile.hpp:316
Definition: file_reader.hpp:6
double priv_flux_decay_factor
profiles decay exponentially to f(sml_outpsi)/decay_factor in priv. flux region
Definition: equil.hpp:94
double p_min
Definition: profile.hpp:31
static constexpr int N_Y
Definition: profile.hpp:174
double xpt2_psi
psi coordinate at 2nd X-point
Definition: equil.hpp:88
CustomLinShape()
Definition: profile.hpp:143
void exit_XGC(std::string msg)
Definition: globals.hpp:37
Shape
Definition: profile.hpp:15
void init_ftn_lin(int num, double *psi, double *var, double psi_min, double dpsi, int n, double *v)
Profile(double val_at_zero_psi, double slope)
Definition: profile.hpp:275
Definition: profile.hpp:24
KOKKOS_INLINE_FUNCTION double value(const MagneticField< DeviceType > &b_field, double psi_in, double r, double z) const
Definition: profile.hpp:319
double xpt_z
z coordinate of 1st X-point
Definition: equil.hpp:83
Definition: profile.hpp:29
Profile(Shape shape_in)
Definition: profile.hpp:286
Profile(const View< double *, CLayout, HostType > &psi, const View< double *, CLayout, HostType > &var)
Definition: profile.hpp:269
double inx[N_X]
Definition: profile.hpp:180
Definition: profile.hpp:16
static constexpr int n
Definition: profile.hpp:30
Kokkos::View< double *, Kokkos::LayoutRight, Device > v
Definition: profile.hpp:34
Kokkos::ViewAllocateWithoutInitializing NoInit
Definition: space_settings.hpp:68
Definition: profile.hpp:171
double p_max
Definition: profile.hpp:32
KOKKOS_INLINE_FUNCTION bool is_in_region_1_or_2(double r, double z, double psi) const
Definition: equil.tpp:82
double priv_flux_decay_width
width for exponential decay in private flux region
Definition: equil.hpp:96
void open(const std::string &input)
Definition: file_reader.hpp:13
double xpt_psi
Psi coordinate of 1st X-point.
Definition: equil.hpp:79