XGC1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
follow_psi_gradients.hpp
Go to the documentation of this file.
1 #ifndef FOLLOW_PSI_GRADIENTS_HPP
2 #define FOLLOW_PSI_GRADIENTS_HPP
3 
4 #include "magnetic_field.hpp"
5 #include "plane.hpp"
6 
7 // RK2 routine to follow the Psi gradient a given distance dPsi in Psi space
8 KOKKOS_INLINE_FUNCTION void get_next_point_gl(const MagneticField<DeviceType>& magnetic_field, double x_in_r, double x_in_z, double psi_in, double dpsi, int direction, double& x_out_r, double& x_out_z, double& dist){
9  constexpr int maxsteps = 20;
10 
11  double x2_r = x_in_r;
12  double x2_z = x_in_z;
13 
14  // Plan for maxsteps/2 iterations -->
15  double dp = direction*dpsi/(maxsteps/2);
16  dist = 0.0;
17 
18  for (int i=0; i<maxsteps; i++){
19  double x1_r = x2_r;
20  double x1_z = x2_z;
21 
22  double psi,dpsi_dr,dpsi_dz;
23 
24  // Predictor step: dl_perp=dpsi/|grad(Psi)|
25  // ==> dl_perp uvector_perp = dpsi*grad(Psi)/|grad(Psi)|^2
26  magnetic_field.psi_bicub.der_one(x1_r, x1_z, psi,dpsi_dr,dpsi_dz );
27  double gradpsi2 = (dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz); // |grad(Psi)|^2
28  x2_r = x1_r + (0.5*dp)*dpsi_dr/gradpsi2;
29  x2_z = x1_z + (0.5*dp)*dpsi_dz/gradpsi2;
30 
31  // I need to implement a fallback in case we are on the separatrix:
32  if(gradpsi2<1.0e-13){
33  DEVICE_PRINTF("get_next_point_gl: Warning, gradient smaller than threshold!\n");
34  }
35 
36  // Corrector step
37  magnetic_field.psi_bicub.der_one(x2_r, x2_z, psi,dpsi_dr,dpsi_dz );
38  gradpsi2 = (dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
39  x2_r = x1_r + dp*dpsi_dr/gradpsi2;
40  x2_z = x1_z + dp*dpsi_dz/gradpsi2;
41 
42  magnetic_field.psi_bicub.der_zero(x2_r, x2_z, psi);
43  if(abs(psi-psi_in)>=dpsi){
44  if(abs(psi-psi_in)>1.5*dpsi){
45  // We have gone too far, go one step back
46  x2_r = x1_r;
47  x2_z = x1_z;
48  } else {
49  dist += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
50  }
51  // exit the loop
52  x_out_r = x2_r;
53  x_out_z = x2_z;
54  break;
55  }
56  dist += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
57  }
58 }
59 
60 // RK2 routine to follow the Psi gradient a given distance dlperp
61 KOKKOS_INLINE_FUNCTION void get_next_point_perp(const MagneticField<DeviceType>& magnetic_field, double x_in_r, double x_in_z, double dlperp, int direction, double& x_out_r, double& x_out_z, double& dist){
62  constexpr int maxsteps = 20;
63 
64  double x2_r = x_in_r;
65  double x2_z = x_in_z;
66 
67  // Plan for maxsteps/2 iterations -->
68  double dlp = direction*dlperp/(maxsteps/2);
69  dist = 0.0;
70  double dlc = 0;
71 
72  for (int i=0; i<maxsteps; i++){
73  double x1_r = x2_r;
74  double x1_z = x2_z;
75 
76  double psi,dpsi_dr,dpsi_dz;
77 
78  // Predictor step
79  magnetic_field.psi_bicub.der_one(x1_r, x1_z, psi,dpsi_dr,dpsi_dz );
80  double gradpsi = sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
81  double inv_gradpsi = 1.0/gradpsi;
82  x2_r = x1_r + (0.5*dlp)*dpsi_dr*inv_gradpsi;
83  x2_z = x1_z + (0.5*dlp)*dpsi_dz*inv_gradpsi;
84 
85  // I need to implement a fallback in case we are on the separatrix:
86  if(gradpsi<1.0e-8 && abs(x1_r - magnetic_field.equil.axis_r)>1.0e-4 && abs(x1_z - magnetic_field.equil.axis_z)>1.0e-9){
87  DEVICE_PRINTF("get_next_point_perp: Warning, gradient smaller than threshold! %1.15e %1.15e %1.15e %1.15e\n",x1_r,x1_z,
88  abs(x1_r - magnetic_field.equil.axis_r), abs(x1_z - magnetic_field.equil.axis_z));
89  }
90 
91  // Corrector step
92  magnetic_field.psi_bicub.der_one(x2_r, x2_z, psi,dpsi_dr,dpsi_dz );
93  gradpsi = sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
94  inv_gradpsi = 1.0/gradpsi;
95  x2_r = x1_r + dlp*dpsi_dr*inv_gradpsi;
96  x2_z = x1_z + dlp*dpsi_dz*inv_gradpsi;
97 
98  dlc += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
99  if(dlc>=dlperp){
100  if(dlc>1.5*dlperp){
101  // We have gone too far, go one step back
102  dist = dlc - sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
103  x2_r = x1_r;
104  x2_z = x1_z;
105  } else {
106  dist = dlc;
107  }
108  // exit the loop
109  x_out_r = x2_r;
110  x_out_z = x2_z;
111  break;
112  }
113  dist = dlc;
114  }
115 }
116 
117 // RK4 routine to follow the flux surface
118 KOKKOS_INLINE_FUNCTION void get_next_point_tang(const MagneticField<DeviceType>& magnetic_field, double x_r_in, double x_z_in, double dltheta, int direction, double& x_r_out, double& x_z_out, double& dist){
119  constexpr int maxsteps = 50;
120 
121  double x_r5 = x_r_in;
122  double x_z5 = x_z_in;
123 
124  // Plan for maxsteps/2 iterations -->
125  double dlt = direction*dltheta/(maxsteps/2);
126  dist = 0.0;
127  double dlc = 0;
128 
129  for (int i=0; i<maxsteps; i++){
130  double inv_gradpsi,psi,dpsi_dr,dpsi_dz;
131 
132  // RK step 1
133  double x_r1 = x_r5;
134  double x_z1 = x_z5;
135  magnetic_field.psi_bicub.der_one(x_r1, x_z1, psi,dpsi_dr,dpsi_dz );
136  inv_gradpsi = 1.0/sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
137  double xp_r1 = -dpsi_dz*inv_gradpsi;
138  double xp_z1 = dpsi_dr*inv_gradpsi;
139 
140  // I need to implement a fallback in case we are on the separatrix:
141  //if (abspolv .lt. 1d-13) then
142  // print *, 'get_next_point_gl: Warning, gradient smaller than threshold!'
143  //endif
144 
145  // RK step 2
146  double x_r2 = x_r1 + (0.5*dlt)*xp_r1;
147  double x_z2 = x_z1 + (0.5*dlt)*xp_z1;
148  magnetic_field.psi_bicub.der_one(x_r2, x_z2, psi,dpsi_dr,dpsi_dz );
149  inv_gradpsi = 1.0/sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
150  double xp_r2 = -dpsi_dz*inv_gradpsi;
151  double xp_z2 = dpsi_dr*inv_gradpsi;
152 
153  // RK step 3
154  double x_r3 = x_r1 + (0.5*dlt)*xp_r2;
155  double x_z3 = x_z1 + (0.5*dlt)*xp_z2;
156  magnetic_field.psi_bicub.der_one(x_r3, x_z3, psi,dpsi_dr,dpsi_dz );
157  inv_gradpsi = 1.0/sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
158  double xp_r3 = -dpsi_dz*inv_gradpsi;
159  double xp_z3 = dpsi_dr*inv_gradpsi;
160 
161  // RK step 4
162  double x_r4 = x_r1 + dlt*xp_r3;
163  double x_z4 = x_z1 + dlt*xp_z3;
164  magnetic_field.psi_bicub.der_one(x_r4, x_z4, psi,dpsi_dr,dpsi_dz );
165  inv_gradpsi = 1.0/sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
166  double xp_r4 = -dpsi_dz*inv_gradpsi;
167  double xp_z4 = dpsi_dr*inv_gradpsi;
168 
169  // Corrector step
170  x_r5 = x_r1 + (dlt/6.0)*(xp_r1 + 2.0*(xp_r2 + xp_r3) + xp_r4);
171  x_z5 = x_z1 + (dlt/6.0)*(xp_z1 + 2.0*(xp_z2 + xp_z3) + xp_z4);
172 
173  dlc += sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
174  if(dlc>=dltheta){
175  if(dlc>1.5*dltheta){
176  // We have gone too far, go one step back
177  dist = dlc - sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
178  x_r5 = x_r1;
179  x_z5 = x_z1;
180  } else {
181  dist = dlc;
182  }
183  // exit the loop
184  x_r_out = x_r5;
185  x_z_out = x_z5;
186  break;
187  }
188  dist = dlc;
189  }
190 }
191 // This routine uses get_next_point_gl to find a valid point in grad(psi)
192 // direction using a psi-difference. It reduces the input distance until it finds a valid point
193 // or the target distance is less than 10% of the input distance
194 KOKKOS_INLINE_FUNCTION void get_node_psi(const MagneticField<DeviceType>& magnetic_field, const Plane<DeviceType>& plane, double x_r, double x_z, double psi_in, double dpsi_in, int direction, double& dist_out, int& itr, double (&p)[3]){
195 
196  double dpsi = dpsi_in;
197  itr = -1;
198 
199  double x2_r = -1.0;
200  double x2_z = -1.0;
201 
202  while(itr<1 && dpsi>=0.1*dpsi_in){
203  get_next_point_gl(magnetic_field, x_r, x_z, psi_in,dpsi,direction, x2_r, x2_z, dist_out);
204  double psi;
205  magnetic_field.psi_bicub.der_zero(x2_r, x2_z, psi);
206  plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
207 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
208  plane.t_coeff_mod(magnetic_field, x2_r,x2_z,psi,itr,p);
209 #endif
210 
211  if(itr<1){
212  dpsi *= 0.4;
213  } else {
214  break;
215  }
216  }
217  if(itr<1){
218  // This is not an error, although it should not happen.
219  // This is likely a node close to the wall.
220  dist_out = 0.0;
221  // if (sml_mype==0) then
222  // print *,'init_grid_deriv: No base point found in Psi-direction'
223  // print '(a11,i3)', 'Direction: ', direction
224  // print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
225  // endif
226  }
227 
228 }
229 
230 // This routine uses get_next_point_perp to find a valid point in grad(psi)
231 // direction using a distance in m. It reduces the input distance until it finds a valid point
232 // or the target distance is less than 10% of the input distance
233 KOKKOS_INLINE_FUNCTION void get_node_perp(const MagneticField<DeviceType>& magnetic_field, const Plane<DeviceType>& plane, double x_r, double x_z, double dlperp_in, int direction, double& dist_out, int& itr, double (&p)[3]){
234 
235  double dlperp = dlperp_in;
236  itr = -1;
237 
238  double x2_r = -1.0;
239  double x2_z = -1.0;
240  while(itr<1 && dlperp>=0.1*dlperp_in){
241  get_next_point_perp(magnetic_field, x_r,x_z,dlperp,direction, x2_r, x2_z, dist_out);
242  double psi;
243  magnetic_field.psi_bicub.der_zero(x2_r, x2_z, psi);
244  plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
245 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
246  plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
247 #endif
248 
249  if(itr<1){
250  dlperp *= 0.4;
251  } else {
252  break;
253  }
254  }
255  if(itr<1){
256  // This is not an error, although it should not happen.
257  // This is likely a node close to the wall.
258  dist_out = 0.0;
259  // if (sml_mype==0) then
260  // print *,'init_grid_deriv: No base point found in Psi-direction'
261  // print '(a11,i3)', 'Direction: ', direction
262  // print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
263  // endif
264  }
265 
266 }
267 
268 // This routine uses get_next_point_tang to find a valid point in grad(theta)
269 // direction using a distance in m. It reduces the input distance until it finds a valid point
270 // or the target distance is less than 10% of the input distance
271 KOKKOS_INLINE_FUNCTION void get_node_theta(const MagneticField<DeviceType>& magnetic_field, const Plane<DeviceType>& plane, double x_r, double x_z, double dltheta_in, int direction, double& dist_out, int& itr, double (&p)[3]){
272  double dltheta = dltheta_in;
273  itr = -1;
274 
275  double x2_r = -1.0;
276  double x2_z = -1.0;
277  while(itr<1 && dltheta>=0.1*dltheta_in){
278  get_next_point_tang(magnetic_field, x_r, x_z, dltheta, direction, x2_r, x2_z, dist_out);
279  double psi;
280  magnetic_field.psi_bicub.der_zero(x2_r, x2_z, psi);
281  plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
282 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
283  plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
284 #endif
285  if(itr<1){
286  dltheta *= 0.4;
287  } else {
288  break;
289  }
290  }
291  if(itr<1){
292  // This is not an error, although it should not happen.
293  // This is likely a node close to the wall.
294  dist_out = 0.0;
295  // if (sml_mype==0) then
296  // print *,'init_grid_deriv: No base point found in theta-direction'
297  // print '(a11,i3)', 'Direction: ', direction
298  // print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
299  // endif
300  }
301 
302 }
303 
304 #endif
subroutine get_node_perp(grid, x, dlperp_in, direction, dist_out, itr, p)
Definition: search.F90:3024
#define DEVICE_PRINTF(...)
Definition: space_settings.hpp:85
Definition: magnetic_field.hpp:12
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:32
KOKKOS_INLINE_FUNCTION void t_coeff_mod(const MagneticField< Device > &magnetic_field, const double r, const double z, const double psiin, const int itr, double(&p)[3]) const
Definition: plane.tpp:139
double axis_z
z coordinate of axis
Definition: equil.hpp:90
KOKKOS_INLINE_FUNCTION void get_node_psi(const MagneticField< DeviceType > &magnetic_field, const Plane< DeviceType > &plane, double x_r, double x_z, double psi_in, double dpsi_in, int direction, double &dist_out, int &itr, double(&p)[3])
Definition: follow_psi_gradients.hpp:194
double axis_r
r coordinate of axis
Definition: equil.hpp:89
KOKKOS_INLINE_FUNCTION void search_tr2_no_precheck(const double r, const double z, int &itr, double(&p)[3]) const
Definition: plane.tpp:13
subroutine get_node_theta(grid, x, dltheta_in, direction, dist_out, itr, p)
Definition: search.F90:3067
Definition: magnetic_field.F90:1
subroutine get_next_point_gl(x_in, psi_in, dpsi, direction, x_out, dist)
Definition: search.F90:3107
subroutine get_next_point_perp(x_in, dlperp, direction, x_out, dist)
Definition: search.F90:3158
Bicub< Device > psi_bicub
The object for interpolating psi (magnetic flux surfaces)
Definition: magnetic_field.hpp:30
subroutine get_next_point_tang(x_in, dltheta, direction, x_out, dist)
Definition: search.F90:3210