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,dpsi_dphi;
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);
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.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;
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,dpsi_dphi;
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;
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.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;
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));
118 double psi,dpsi_dr,dpsi_dz,dpsi_dphi;
120 double inv_gradpsi = 1.0/sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
122 rp = -dpsi_dz*inv_gradpsi;
123 zp = dpsi_dr*inv_gradpsi;
133 constexpr
int maxsteps = 50;
135 double x_r5 = x_r_in;
136 double x_z5 = x_z_in;
139 double dlt = direction*dltheta/(maxsteps/2);
143 for (
int i=0; i<maxsteps; i++){
151 double x_r2 = x_r1 + (0.5*dlt)*xp_r1;
152 double x_z2 = x_z1 + (0.5*dlt)*xp_z1;
157 double x_r3 = x_r1 + (0.5*dlt)*xp_r2;
158 double x_z3 = x_z1 + (0.5*dlt)*xp_z2;
163 double x_r4 = x_r1 + dlt*xp_r3;
164 double x_z4 = x_z1 + dlt*xp_z3;
169 x_r5 = x_r1 + (dlt/6.0)*(xp_r1 + 2.0*(xp_r2 + xp_r3) + xp_r4);
170 x_z5 = x_z1 + (dlt/6.0)*(xp_z1 + 2.0*(xp_z2 + xp_z3) + xp_z4);
172 dlc += sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
176 dist = dlc - sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
195 double dpsi = dpsi_in;
201 while(itr<1 && dpsi>=0.1*dpsi_in){
205 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
233 double dlperp = dlperp_in;
238 while(itr<1 && dlperp>=0.1*dlperp_in){
242 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
269 double dltheta = dltheta_in;
274 while(itr<1 && dltheta>=0.1*dltheta_in){
278 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
299 double dltheta = dltheta_in;
308 while(itr<1 && dltheta>=0.1*dltheta_in){
312 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
Definition: magnetic_field.hpp:12
KOKKOS_INLINE_FUNCTION void search_tr2_no_precheck(const double r, const double z, int &itr, double(&p)[3]) const
Definition: plane.tpp:13
double plane_phi
phi coordinate of this plane
Definition: plane.hpp:261
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
KOKKOS_INLINE_FUNCTION void get_rz_coordinates(const int inode, double &r, double &z) const
Definition: plane.tpp:1003
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:193
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)
Definition: follow_psi_gradients.hpp:61
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)
Definition: follow_psi_gradients.hpp:8
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])
Definition: follow_psi_gradients.hpp:268
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])
Definition: follow_psi_gradients.hpp:231
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)
Definition: follow_psi_gradients.hpp:132
KOKKOS_INLINE_FUNCTION void single_tang_step(const MagneticField< DeviceType > &magnetic_field, double r, double z, double phi, double &rp, double &zp)
Definition: follow_psi_gradients.hpp:117
integer, parameter dp
Definition: module.F90:2453
real(8), parameter phi
Definition: load_balance_constraint_mod.F90:18
Definition: magnetic_field.F90:1
#define DEVICE_PRINTF(...)
Definition: space_settings.hpp:86