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 phi, 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,dpsi_dphi;
23 
24  // Predictor step: dl_perp=dpsi/|grad(Psi)|
25  // ==> dl_perp uvector_perp = dpsi*grad(Psi)/|grad(Psi)|^2
26  magnetic_field.get_psi_and_derivs(x1_r, x1_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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.get_psi_and_derivs(x2_r, x2_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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  psi = magnetic_field.get_psi(x2_r, x2_z, phi);
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 phi, 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,dpsi_dphi;
77 
78  // Predictor step
79  magnetic_field.get_psi_and_derivs(x1_r, x1_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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.get_psi_and_derivs(x2_r, x2_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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 phi, 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,dpsi_dphi;
131 
132  // RK step 1
133  double x_r1 = x_r5;
134  double x_z1 = x_z5;
135  magnetic_field.get_psi_and_derivs(x_r1, x_z1, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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.get_psi_and_derivs(x_r2, x_z2, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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.get_psi_and_derivs(x_r3, x_z3, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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.get_psi_and_derivs(x_r4, x_z4, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
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 phi, 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, phi, psi_in,dpsi,direction, x2_r, x2_z, dist_out);
204  double psi = magnetic_field.get_psi(x2_r, x2_z, phi);
205  plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
206 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
207  plane.t_coeff_mod(magnetic_field, x2_r,x2_z,psi,itr,p);
208 #endif
209 
210  if(itr<1){
211  dpsi *= 0.4;
212  } else {
213  break;
214  }
215  }
216  if(itr<1){
217  // This is not an error, although it should not happen.
218  // This is likely a node close to the wall.
219  dist_out = 0.0;
220  // if (sml_mype==0) then
221  // print *,'init_grid_deriv: No base point found in Psi-direction'
222  // print '(a11,i3)', 'Direction: ', direction
223  // print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
224  // endif
225  }
226 
227 }
228 
229 // This routine uses get_next_point_perp to find a valid point in grad(psi)
230 // direction using a distance in m. It reduces the input distance until it finds a valid point
231 // or the target distance is less than 10% of the input distance
232 KOKKOS_INLINE_FUNCTION void get_node_perp(const MagneticField<DeviceType>& magnetic_field, const Plane<DeviceType>& plane, double x_r, double x_z, double phi, double dlperp_in, int direction, double& dist_out, int& itr, double (&p)[3]){
233 
234  double dlperp = dlperp_in;
235  itr = -1;
236 
237  double x2_r = -1.0;
238  double x2_z = -1.0;
239  while(itr<1 && dlperp>=0.1*dlperp_in){
240  get_next_point_perp(magnetic_field, x_r,x_z,phi, dlperp,direction, x2_r, x2_z, dist_out);
241  double psi = magnetic_field.get_psi(x2_r, x2_z, phi);
242  plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
243 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
244  plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
245 #endif
246 
247  if(itr<1){
248  dlperp *= 0.4;
249  } else {
250  break;
251  }
252  }
253  if(itr<1){
254  // This is not an error, although it should not happen.
255  // This is likely a node close to the wall.
256  dist_out = 0.0;
257  // if (sml_mype==0) then
258  // print *,'init_grid_deriv: No base point found in Psi-direction'
259  // print '(a11,i3)', 'Direction: ', direction
260  // print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
261  // endif
262  }
263 
264 }
265 
266 // This routine uses get_next_point_tang to find a valid point in grad(theta)
267 // direction using a distance in m. It reduces the input distance until it finds a valid point
268 // or the target distance is less than 10% of the input distance
269 KOKKOS_INLINE_FUNCTION void get_node_theta(const MagneticField<DeviceType>& magnetic_field, const Plane<DeviceType>& plane, double x_r, double x_z, double phi, double dltheta_in, int direction, double& dist_out, int& itr, double (&p)[3]){
270  double dltheta = dltheta_in;
271  itr = -1;
272 
273  double x2_r = -1.0;
274  double x2_z = -1.0;
275  while(itr<1 && dltheta>=0.1*dltheta_in){
276  get_next_point_tang(magnetic_field, x_r, x_z, phi, dltheta, direction, x2_r, x2_z, dist_out);
277  double psi = magnetic_field.get_psi(x2_r, x2_z, phi);
278  plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
279 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
280  plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
281 #endif
282  if(itr<1){
283  dltheta *= 0.4;
284  } else {
285  break;
286  }
287  }
288  if(itr<1){
289  // This is not an error, although it should not happen.
290  // This is likely a node close to the wall.
291  dist_out = 0.0;
292  // if (sml_mype==0) then
293  // print *,'init_grid_deriv: No base point found in theta-direction'
294  // print '(a11,i3)', 'Direction: ', direction
295  // print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
296  // endif
297  }
298 
299 }
300 
301 #endif
KOKKOS_INLINE_FUNCTION void get_psi_and_derivs(double r, double z, double phi, double &psi, double &dpsidr, double &dpsidz, double &dpsidphi) const
Definition: magnetic_field.tpp:121
KOKKOS_INLINE_FUNCTION void get_psi(const SimdVector &v, Simd< double > &psi_out) const
Definition: magnetic_field.tpp:105
subroutine get_node_perp(grid, x, dlperp_in, direction, dist_out, itr, p)
Definition: search.F90:2964
#define DEVICE_PRINTF(...)
Definition: space_settings.hpp:85
Definition: magnetic_field.hpp:13
Equilibrium equil
The object containing information about the magnetic equilibrium.
Definition: magnetic_field.hpp:45
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 z
Definition: grid_structs.hpp:30
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:3007
double r
Definition: grid_structs.hpp:29
Definition: magnetic_field.F90:1
subroutine get_next_point_gl(x_in, psi_in, dpsi, direction, x_out, dist)
Definition: search.F90:3047
subroutine get_next_point_perp(x_in, dlperp, direction, x_out, dist)
Definition: search.F90:3098
RZPair axis
Definition: equil.hpp:86
subroutine get_next_point_tang(x_in, dltheta, direction, x_out, dist)
Definition: search.F90:3150
KOKKOS_INLINE_FUNCTION void get_node_psi(const MagneticField< DeviceType > &magnetic_field, const Plane< DeviceType > &plane, double x_r, double x_z, double phi, double psi_in, double dpsi_in, int direction, double &dist_out, int &itr, double(&p)[3])
Definition: follow_psi_gradients.hpp:194