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;
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");
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 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){
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;
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,
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,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;
146 double x_r2 = x_r1 + (0.5*dlt)*xp_r1;
147 double x_z2 = x_z1 + (0.5*dlt)*xp_z1;
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;
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;
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));
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]){
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, phi, psi_in,dpsi,direction, x2_r, x2_z, dist_out);
204 double psi = magnetic_field.
get_psi(x2_r, x2_z, phi);
206 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
207 plane.
t_coeff_mod(magnetic_field, x2_r,x2_z,psi,itr,p);
234 double dlperp = dlperp_in;
239 while(itr<1 && dlperp>=0.1*dlperp_in){
241 double psi = magnetic_field.
get_psi(x2_r, x2_z, phi);
243 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
244 plane.
t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
270 double dltheta = dltheta_in;
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);
279 #if defined(XGCA) || defined(NEOCLASSICAL_TEST)
280 plane.
t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
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