follow_psi_gradients.hpp

Go to the documentation of this file.

#ifndef FOLLOW_PSI_GRADIENTS_HPP
#define FOLLOW_PSI_GRADIENTS_HPP
 
#include "magnetic_field.hpp"
#include "plane.hpp"
 
// RK2 routine to follow the Psi gradient a given distance dPsi in Psi space
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){
    constexpr int maxsteps = 20;
 
    double x2_r = x_in_r;
    double x2_z = x_in_z;
 
    // Plan for maxsteps/2 iterations -->
    double dp = direction*dpsi/(maxsteps/2);
    dist = 0.0;
 
    for (int i=0; i<maxsteps; i++){
        double x1_r = x2_r;
        double x1_z = x2_z;
 
        double psi,dpsi_dr,dpsi_dz,dpsi_dphi;
 
        // Predictor step: dl_perp=dpsi/|grad(Psi)|
        // ==> dl_perp uvector_perp = dpsi*grad(Psi)/|grad(Psi)|^2
        magnetic_field.get_psi_and_derivs(x1_r, x1_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
        double gradpsi2 = (dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz); // |grad(Psi)|^2
        x2_r = x1_r + (0.5*dp)*dpsi_dr/gradpsi2;
        x2_z = x1_z + (0.5*dp)*dpsi_dz/gradpsi2;
 
        // I need to implement a fallback in case we are on the separatrix:
        if(gradpsi2<1.0e-13){
            DEVICE_PRINTF("get_next_point_gl: Warning, gradient smaller than threshold!\n");
        }
 
        // Corrector step
        magnetic_field.get_psi_and_derivs(x2_r, x2_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
        gradpsi2 = (dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
        x2_r = x1_r + dp*dpsi_dr/gradpsi2;
        x2_z = x1_z + dp*dpsi_dz/gradpsi2;
 
        psi = magnetic_field.get_psi(x2_r, x2_z, phi);
        if(abs(psi-psi_in)>=dpsi){
            if(abs(psi-psi_in)>1.5*dpsi){
                // We have gone too far, go one step back
                x2_r = x1_r;
                x2_z = x1_z;
            } else {
                dist += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
            }
            // exit the loop
            x_out_r = x2_r;
            x_out_z = x2_z;
            break;
        }
        dist += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
    }
}
 
// RK2 routine to follow the Psi gradient a given distance dlperp
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){
    constexpr int maxsteps = 20;
 
    double x2_r = x_in_r;
    double x2_z = x_in_z;
 
    // Plan for maxsteps/2 iterations -->
    double dlp = direction*dlperp/(maxsteps/2);
    dist = 0.0;
    double dlc = 0;
 
    for (int i=0; i<maxsteps; i++){
        double x1_r = x2_r;
        double x1_z = x2_z;
 
        double psi,dpsi_dr,dpsi_dz,dpsi_dphi;
 
        // Predictor step
        magnetic_field.get_psi_and_derivs(x1_r, x1_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
        double gradpsi = sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
        double inv_gradpsi = 1.0/gradpsi;
        x2_r = x1_r + (0.5*dlp)*dpsi_dr*inv_gradpsi;
        x2_z = x1_z + (0.5*dlp)*dpsi_dz*inv_gradpsi;
 
        // I need to implement a fallback in case we are on the separatrix:
        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){
            DEVICE_PRINTF("get_next_point_perp: Warning, gradient smaller than threshold! %1.15e %1.15e %1.15e %1.15e\n",x1_r,x1_z,
                          abs(x1_r - magnetic_field.equil.axis.r), abs(x1_z - magnetic_field.equil.axis.z));
        }
 
        // Corrector step
        magnetic_field.get_psi_and_derivs(x2_r, x2_z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
        gradpsi = sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
        inv_gradpsi = 1.0/gradpsi;
        x2_r = x1_r + dlp*dpsi_dr*inv_gradpsi;
        x2_z = x1_z + dlp*dpsi_dz*inv_gradpsi;
 
        dlc += sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
        if(dlc>=dlperp){
            if(dlc>1.5*dlperp){
                // We have gone too far, go one step back
                dist = dlc - sqrt((x2_r-x1_r)*(x2_r-x1_r) + (x2_z-x1_z)*(x2_z-x1_z));
                x2_r = x1_r;
                x2_z = x1_z;
            } else {
                dist = dlc;
            }
            // exit the loop
            x_out_r = x2_r;
            x_out_z = x2_z;
            break;
        }
        dist = dlc;
    }
}
 
KOKKOS_INLINE_FUNCTION void single_tang_step(const MagneticField<DeviceType>& magnetic_field, double r, double z, double phi, double& rp, double& zp){
    double psi,dpsi_dr,dpsi_dz,dpsi_dphi;
    magnetic_field.get_psi_and_derivs(r, z, phi, psi,dpsi_dr,dpsi_dz,dpsi_dphi );
    double inv_gradpsi = 1.0/sqrt(dpsi_dr*dpsi_dr + dpsi_dz*dpsi_dz);
    //double inv_gradpsi = 1.0/gradpsi;
    rp = -dpsi_dz*inv_gradpsi;
    zp = dpsi_dr*inv_gradpsi;
 
    // If gradpsi is small, set to derivatives to 0
    //constexpr double GRADPSI_MIN = 1.0e-12;
    //rp = gradpsi<GRADPSI_MIN ? 0.0 : rp;
    //zp = gradpsi<GRADPSI_MIN ? 0.0 : zp;
}
 
// RK4 routine to follow the flux surface
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){
    constexpr int maxsteps = 50;
 
    double x_r5 = x_r_in;
    double x_z5 = x_z_in;
 
    // Plan for maxsteps/2 iterations -->
    double dlt = direction*dltheta/(maxsteps/2);
    dist = 0.0;
    double dlc = 0;
 
    for (int i=0; i<maxsteps; i++){
        // RK step 1
        double x_r1 = x_r5;
        double x_z1 = x_z5;
        double xp_r1, xp_z1;
        single_tang_step(magnetic_field, x_r1, x_z1, phi, xp_r1, xp_z1);
 
        // RK step 2
        double x_r2 = x_r1 + (0.5*dlt)*xp_r1;
        double x_z2 = x_z1 + (0.5*dlt)*xp_z1;
        double xp_r2, xp_z2;
        single_tang_step(magnetic_field, x_r2, x_z2, phi, xp_r2, xp_z2);
 
        // RK step 3
        double x_r3 = x_r1 + (0.5*dlt)*xp_r2;
        double x_z3 = x_z1 + (0.5*dlt)*xp_z2;
        double xp_r3, xp_z3;
        single_tang_step(magnetic_field, x_r3, x_z3, phi, xp_r3, xp_z3);
 
        // RK step 4
        double x_r4 = x_r1 + dlt*xp_r3;
        double x_z4 = x_z1 + dlt*xp_z3;
        double xp_r4, xp_z4;
        single_tang_step(magnetic_field, x_r4, x_z4, phi, xp_r4, xp_z4);
 
        // Corrector step
        x_r5 = x_r1 + (dlt/6.0)*(xp_r1 + 2.0*(xp_r2 + xp_r3) + xp_r4);
        x_z5 = x_z1 + (dlt/6.0)*(xp_z1 + 2.0*(xp_z2 + xp_z3) + xp_z4);
 
        dlc += sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
        if(dlc>=dltheta){
            if(dlc>1.5*dltheta){
                // We have gone too far, go one step back
                dist = dlc - sqrt((x_r5-x_r1)*(x_r5-x_r1) + (x_z5-x_z1)*(x_z5-x_z1));
                x_r5 = x_r1;
                x_z5 = x_z1;
            } else {
                dist = dlc;
            }
            // exit the loop
            x_r_out = x_r5;
            x_z_out = x_z5;
            break;
        }
        dist = dlc;
    }
}
// This routine uses get_next_point_gl to find a valid point in grad(psi)
// direction using a psi-difference. It reduces the input distance until it finds a valid point
// or the target distance is less than 10% of the input distance
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]){
 
    double dpsi = dpsi_in;
    itr = -1;
 
    double x2_r = -1.0;
    double x2_z = -1.0;
 
    while(itr<1 && dpsi>=0.1*dpsi_in){
        get_next_point_gl(magnetic_field, x_r, x_z, plane.plane_phi, psi_in,dpsi,direction, x2_r, x2_z, dist_out);
        double psi = magnetic_field.get_psi(x2_r, x2_z, plane.plane_phi);
        plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
#if defined(XGCA) || defined(NEOCLASSICAL_TEST)
        plane.t_coeff_mod(magnetic_field, x2_r,x2_z,psi,itr,p);
#endif
 
        if(itr<1){
            dpsi *= 0.4;
        } else {
            break;
        }
    }
    if(itr<1){
        // This is not an error, although it should not happen.
        // This is likely a node close to the wall.
        dist_out = 0.0;
        //    if (sml_mype==0) then
        //      print *,'init_grid_deriv: No base point found in Psi-direction'
        //      print '(a11,i3)', 'Direction: ', direction
        //      print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
        //    endif
    }
 
}
 
// This routine uses get_next_point_perp to find a valid point in grad(psi)
// direction using a distance in m. It reduces the input distance until it finds a valid point
// or the target distance is less than 10% of the input distance
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]){
 
    double dlperp = dlperp_in;
    itr = -1;
 
    double x2_r = -1.0;
    double x2_z = -1.0;
    while(itr<1 && dlperp>=0.1*dlperp_in){
        get_next_point_perp(magnetic_field, x_r,x_z,plane.plane_phi, dlperp,direction, x2_r, x2_z, dist_out);
        double psi = magnetic_field.get_psi(x2_r, x2_z, plane.plane_phi);
        plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
#if defined(XGCA) || defined(NEOCLASSICAL_TEST)
        plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
#endif
 
        if(itr<1){
            dlperp *= 0.4;
        } else {
            break;
        }
    }
    if(itr<1){
        // This is not an error, although it should not happen.
        // This is likely a node close to the wall.
        dist_out = 0.0;
        //    if (sml_mype==0) then
        //      print *,'init_grid_deriv: No base point found in Psi-direction'
        //      print '(a11,i3)', 'Direction: ', direction
        //      print '(a13,2f16.6)', 'Coordinates: ',x(1),x(2)
        //    endif
    }
 
}
 
// This routine uses get_next_point_tang to find a valid point in grad(theta)
// direction using a distance in m. It reduces the input distance until it finds a valid point
// or the target distance is less than 10% of the input distance
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]){
    double dltheta = dltheta_in;
    itr = -1;
 
    double x2_r = -1.0;
    double x2_z = -1.0;
    while(itr<1 && dltheta>=0.1*dltheta_in){
        get_next_point_tang(magnetic_field, x_r, x_z, plane.plane_phi, dltheta, direction, x2_r, x2_z, dist_out);
        double psi = magnetic_field.get_psi(x2_r, x2_z, plane.plane_phi);
        plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
#if defined(XGCA) || defined(NEOCLASSICAL_TEST)
        plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
#endif
        if(itr<1){
            dltheta *= 0.4;
        } else {
            break;
        }
    }
    if(itr<1){
        // This is not an error, although it should not happen.
        // This is likely a node close to the wall.
        dist_out = 0.0;
    }
}
 
// This routine uses get_next_point_tang to find a valid point in grad(theta)
// direction using a distance in m. It reduces the input distance until it finds a valid point
// or the target distance is less than 10% of the input distance
// Same as above, but input coordinates are provided by itr0 and p0 rather than x_r and x_z
KOKKOS_INLINE_FUNCTION void get_node_theta(const MagneticField<DeviceType>& magnetic_field, const Plane<DeviceType>& plane, int& itr0, double (&p0)[3], double dltheta_in, int direction, double& dist_out, int& itr, double (&p)[3]){
    double dltheta = dltheta_in;
    itr = -1;
    
    // Get rz coordinates since itr0 and p0 are provided rather than x_r and x_z directly
    double x_r, x_z;
    plane.get_rz_coordinates(itr0, p0, x_r, x_z);
 
    double x2_r = -1.0;
    double x2_z = -1.0;
    while(itr<1 && dltheta>=0.1*dltheta_in){
        get_next_point_tang(magnetic_field, x_r, x_z, plane.plane_phi, dltheta, direction, x2_r, x2_z, dist_out);
        double psi = magnetic_field.get_psi(x2_r, x2_z, plane.plane_phi);
        plane.search_tr2_no_precheck(x2_r,x2_z,itr,p);
#if defined(XGCA) || defined(NEOCLASSICAL_TEST)
        plane.t_coeff_mod(magnetic_field,x2_r,x2_z,psi,itr,p);
#endif  
        if(itr<1){
            dltheta *= 0.5;
        } else {
            break;
        }
    }
    if(itr<1){
        // This is not an error, although it should not happen.
        // This is likely a node close to the wall.
        dist_out = 0.0;
    }
}
 
#endif