1 #ifndef FOLLOW_PSI_GRADIENTS_HPP
2 #define FOLLOW_PSI_GRADIENTS_HPP
9 constexpr
int maxsteps = 20;
15 double dp = direction*dpsi/(maxsteps/2);
18 for (
int i=0; i<maxsteps; i++){
22 double psi,dpsi_dr,dpsi_dz;
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);
28 x2_r = x1_r + (0.5*dp)*dpsi_dr/gradpsi2;
29 x2_z = x1_z + (0.5*dp)*dpsi_dz/gradpsi2;
33 DEVICE_PRINTF(
"get_next_point_gl: Warning, gradient smaller than threshold!\n");
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;
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){
49 dist += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
56 dist += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
62 constexpr
int maxsteps = 20;
68 double dlp = direction*dlperp/(maxsteps/2);
72 for (
int i=0; i<maxsteps; i++){
76 double psi,dpsi_dr,dpsi_dz;
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;
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,
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;
98 dlc += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
102 dist = dlc - sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
119 constexpr
int maxsteps = 50;
121 double x_r5 = x_r_in;
122 double x_z5 = x_z_in;
125 double dlt = direction*dltheta/(maxsteps/2);
129 for (
int i=0; i<maxsteps; i++){
130 double inv_gradpsi,psi,dpsi_dr,dpsi_dz;
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;
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;
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;
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;
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);
173 dlc += sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
177 dist = dlc - sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
196 double dpsi = dpsi_in;
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);
205 magnetic_field.
psi_bicub.der_zero(x2_r, x2_z, psi);
207 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
208 grid.
t_coeff_mod(magnetic_field, x2_r,x2_z,psi,itr,p);
235 double dlperp = dlperp_in;
240 while(itr<1 && dlperp>=0.1*dlperp_in){
243 magnetic_field.
psi_bicub.der_zero(x2_r, x2_z, psi);
245 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
246 grid.
t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
272 double dltheta = dltheta_in;
277 while(itr<1 && dltheta>=0.1*dltheta_in){
280 magnetic_field.
psi_bicub.der_zero(x2_r, x2_z, psi);
282 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
283 grid.
t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
subroutine get_node_perp(grid, x, dlperp_in, direction, dist_out, itr, p)
Definition: search.F90:3512
#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: grid.tpp:359
double axis_z
z coordinate of axis
Definition: equil.hpp:95
double axis_r
r coordinate of axis
Definition: equil.hpp:94
KOKKOS_INLINE_FUNCTION void search_tr2_no_precheck(const double r, const double z, int &itr, double(&p)[3]) const
Definition: grid.tpp:233
subroutine get_node_theta(grid, x, dltheta_in, direction, dist_out, itr, p)
Definition: search.F90:3555
KOKKOS_INLINE_FUNCTION void get_node_psi(const MagneticField< DeviceType > &magnetic_field, const Grid< DeviceType > &grid, 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
Definition: magnetic_field.F90:1
subroutine get_next_point_gl(x_in, psi_in, dpsi, direction, x_out, dist)
Definition: search.F90:3595
subroutine get_next_point_perp(x_in, dlperp, direction, x_out, dist)
Definition: search.F90:3646
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:3698