Fieldline following and rotational transform calculation available in 'develop'

[zhucaoxiang]

Hi, I have pushed commits in the branch 'develop'. Now FOCUS can follow fieldlines and calculate rotational transform for vacuum magnetic field produced by coils. Here are some instructions you might need.

1. How to use?
Set case_postproc = 3 and provide appropriate input arguments. Here is the list of parameters.

 pp_phi         =  0.000D+00		! toroidal plane for poincare plots, cylindrical angle phi = pp_phi*Pi
 pp_raxis       =  0.000D+00       	! pp_raxis, pp_zaxis are initial guesses for magnetic axis at the specified toroidal angle
 pp_zaxis       =  0.000D+00            ! If both zero, FOCUS will take the geometric center as initial guess
 pp_rmax        =  0.000D+00		! pp_rmax, pp_zmax are the upper bounds for performing fieldline tracing
 pp_zmax        =  0.000D+00		! FOCUS will start follow fieldlines at interpolation between (pp_raxis, pp_zaxis) and (pp_rmax, pp_zmax)
 pp_ns          =  10			! number of following fieldlines
 pp_maxiter     =  1000			! number of periods for each fieldline following
 pp_xtol        =  1.000D-06		! tolarence of ODE solver during fieldline fowllowing

By default, FOCUS will find the magnetic axis with an initial guess of [Ra,Za]=[R(0,phi) + R(Pi,phi))/2, R(0,phi) + R(Pi,phi))/2 ] at phi=0.0. Then interpolate pp_ns=10 points between [Ra, Za] and [Rmax, Zmax] (if both zero, will set to [R(0,phi), Z(0,phi)]). Afterwards, it will following the fieldlines for 'pp_maxiter' times and calculate the rotational transform.

The following is piece of screening output with settings of pp_ns = 50, pp_phi = 0.5, pp_rmax = 1.4.

 --------------------------------------------
findaxis: Finding axis at phi =   1.57 with (R,Z) = (  8.33002E-01,-2.57496E-08 ).
findaxis: info=1, relative error between two consecutive iterates is at most xtol.
poinplot: following fieldlines between ( 8.33002E-01,-2.57496E-08 ) and ( 1.40000E+00, 0.00000E+00 )
        : order=     1 ; myid=     1 ; (R,Z)=( 8.44342E-01,-2.52346E-08 ) ; iota= 3.88786E-01 ; niter=  1000 .
        : order=     2 ; myid=     2 ; (R,Z)=( 8.55682E-01,-2.47196E-08 ) ; iota= 3.88129E-01 ; niter=  1000 .
        : order=     7 ; myid=     7 ; (R,Z)=( 9.12381E-01,-2.21447E-08 ) ; iota= 3.81505E-01 ; niter=  1000 .
        : order=     3 ; myid=     3 ; (R,Z)=( 8.67022E-01,-2.42046E-08 ) ; iota= 3.87273E-01 ; niter=  1000 .
        : order=     5 ; myid=     5 ; (R,Z)=( 8.89702E-01,-2.31746E-08 ) ; iota= 3.84855E-01 ; niter=  1000 .
        : order=     4 ; myid=     4 ; (R,Z)=( 8.78362E-01,-2.36896E-08 ) ; iota= 3.86127E-01 ; niter=  1000 .
        : order=     6 ; myid=     6 ; (R,Z)=( 9.01041E-01,-2.26597E-08 ) ; iota= 3.83171E-01 ; niter=  1000 .
        : order=     8 ; myid=     0 ; (R,Z)=( 9.23721E-01,-2.16297E-08 ) ; iota= 3.79734E-01 ; niter=  1000 .
        : order=     9 ; myid=     1 ; (R,Z)=( 9.35061E-01,-2.11147E-08 ) ; iota= 3.77851E-01 ; niter=  1000 .
        : order=    10 ; myid=     2 ; (R,Z)=( 9.46401E-01,-2.05997E-08 ) ; iota= 3.75916E-01 ; niter=  1000 .
        : order=    11 ; myid=     3 ; (R,Z)=( 9.57741E-01,-2.00847E-08 ) ; iota= 3.73972E-01 ; niter=  1000 .
        : order=    12 ; myid=     4 ; (R,Z)=( 9.69081E-01,-1.95697E-08 ) ; iota= 3.72034E-01 ; niter=  1000 .
        : order=    13 ; myid=     5 ; (R,Z)=( 9.80421E-01,-1.90547E-08 ) ; iota= 3.70142E-01 ; niter=  1000 .
        : order=    15 ; myid=     7 ; (R,Z)=( 1.00310E+00,-1.80247E-08 ) ; iota= 3.66726E-01 ; niter=  1000 .
        : order=    14 ; myid=     6 ; (R,Z)=( 9.91761E-01,-1.85397E-08 ) ; iota= 3.68386E-01 ; niter=  1000 .
        : order=    16 ; myid=     0 ; (R,Z)=( 1.01444E+00,-1.75097E-08 ) ; iota= 3.64995E-01 ; niter=  1000 .
        : order=    18 ; myid=     2 ; (R,Z)=( 1.03712E+00,-1.64798E-08 ) ; iota= 3.62079E-01 ; niter=  1000 .
        : order=    17 ; myid=     1 ; (R,Z)=( 1.02578E+00,-1.69947E-08 ) ; iota= 3.63685E-01 ; niter=  1000 .
        : order=    19 ; myid=     3 ; (R,Z)=( 1.04846E+00,-1.59648E-08 ) ; iota= 3.60787E-01 ; niter=  1000 .
        : order=    20 ; myid=     4 ; (R,Z)=( 1.05980E+00,-1.54498E-08 ) ; iota= 3.59449E-01 ; niter=  1000 .
        : order=    21 ; myid=     5 ; (R,Z)=( 1.07114E+00,-1.49348E-08 ) ; iota= 3.58292E-01 ; niter=  1000 .
        : order=    22 ; myid=     6 ; (R,Z)=( 1.08248E+00,-1.44198E-08 ) ; iota= 3.57252E-01 ; niter=  1000 .
        : order=    23 ; myid=     7 ; (R,Z)=( 1.09382E+00,-1.39048E-08 ) ; iota= 3.56341E-01 ; niter=  1000 .
        : order=    24 ; myid=     0 ; (R,Z)=( 1.10516E+00,-1.33898E-08 ) ; iota= 3.55539E-01 ; niter=  1000 .
        : order=    26 ; myid=     2 ; (R,Z)=( 1.12784E+00,-1.23598E-08 ) ; iota= 3.54129E-01 ; niter=  1000 .
        : order=    25 ; myid=     1 ; (R,Z)=( 1.11650E+00,-1.28748E-08 ) ; iota= 3.54846E-01 ; niter=  1000 .
        : order=    27 ; myid=     3 ; (R,Z)=( 1.13918E+00,-1.18448E-08 ) ; iota= 3.53689E-01 ; niter=  1000 .
        : order=    28 ; myid=     4 ; (R,Z)=( 1.15052E+00,-1.13298E-08 ) ; iota= 3.53187E-01 ; niter=  1000 .
        : order=    29 ; myid=     5 ; (R,Z)=( 1.16186E+00,-1.08148E-08 ) ; iota= 3.52688E-01 ; niter=  1000 .
        : order=    30 ; myid=     6 ; (R,Z)=( 1.17320E+00,-1.02998E-08 ) ; iota= 3.52278E-01 ; niter=  1000 .
        : order=    31 ; myid=     7 ; (R,Z)=( 1.18454E+00,-9.78485E-09 ) ; iota= 3.52014E-01 ; niter=  1000 .
        : order=    32 ; myid=     0 ; (R,Z)=( 1.19588E+00,-9.26986E-09 ) ; iota= 3.51826E-01 ; niter=  1000 .
        : order=    34 ; myid=     2 ; (R,Z)=( 1.21856E+00,-8.23988E-09 ) ; iota= 3.51271E-01 ; niter=  1000 .
        : order=    33 ; myid=     1 ; (R,Z)=( 1.20722E+00,-8.75487E-09 ) ; iota= 3.51518E-01 ; niter=  1000 .
        : order=    35 ; myid=     3 ; (R,Z)=( 1.22990E+00,-7.72488E-09 ) ; iota= 3.51111E-01 ; niter=  1000 .
        : order=    36 ; myid=     4 ; (R,Z)=( 1.24124E+00,-7.20989E-09 ) ; iota= 3.50933E-01 ; niter=  1000 .
        : order=    37 ; myid=     5 ; (R,Z)=( 1.25258E+00,-6.69490E-09 ) ; iota= 3.50759E-01 ; niter=  1000 .
        : order=    38 ; myid=     6 ; (R,Z)=( 1.26392E+00,-6.17991E-09 ) ; iota= 3.50525E-01 ; niter=  1000 .
        : order=    39 ; myid=     7 ; (R,Z)=( 1.27526E+00,-5.66491E-09 ) ; iota= 3.50189E-01 ; niter=  1000 .
        : order=    40 ; myid=     0 ; (R,Z)=( 1.28660E+00,-5.14992E-09 ) ; iota= 3.49959E-01 ; niter=  1000 .
        : order=    43 ; myid=     3 ; (R,Z)=( 1.32062E+00,-3.60495E-09 ) ; iota= 1.51195E-01 ; niter=  1000 .
        : order=    41 ; myid=     1 ; (R,Z)=( 1.29794E+00,-4.63493E-09 ) ; iota= 3.49558E-01 ; niter=  1000 .
        : order=    42 ; myid=     2 ; (R,Z)=( 1.30928E+00,-4.11994E-09 ) ; iota= 3.49067E-01 ; niter=  1000 .
        : order=    47 ; myid=     7 ; (R,Z)=( 1.36598E+00,-1.54498E-09 ) ; iota= 3.47503E-01 ; niter=  1000 .
        : order=    44 ; myid=     4 ; (R,Z)=( 1.33196E+00,-3.08995E-09 ) ; iota= 3.47897E-01 ; niter=  1000 .
        : order=    45 ; myid=     5 ; (R,Z)=( 1.34330E+00,-2.57496E-09 ) ; iota= 3.47884E-01 ; niter=  1000 .
        : order=    46 ; myid=     6 ; (R,Z)=( 1.35464E+00,-2.05997E-09 ) ; iota= 3.47860E-01 ; niter=  1000 .
        : order=    49 ; myid=     1 ; (R,Z)=( 1.38866E+00,-5.14992E-10 ) ; iota= 2.71804E-02 ; niter=  1000 .
        : order=    50 ; myid=     2 ; (R,Z)=( 1.40000E+00, 0.00000E+00 ) ; iota= 5.19176E-03 ; niter=  1000 .
        : order=    48 ; myid=     0 ; (R,Z)=( 1.37732E+00,-1.02998E-09 ) ; iota= 3.36773E-01 ; niter=    75 .
focus   : Post-processing took      1 M     27 S.
 -------------------------------------------------------------

2. Data processing
The calculated data are store in hdf5 file with the following variables.

ppr(1:pp_ns, 0:pp_maxiter)
ppz(1:pp_ns, 0:pp_maxiter)   
iota(1:pp_ns)

Here is an example of plotting scripts using python (in my personal package coilpy).

def plot_focus_poincare(focus_data, color='k', dotsize=0.1):
    '''
    Poincare plots in FOCUS
    '''
    # poincare plots
    for i in range(focus_data.pp_ns):
        plt.scatter(focus_data.ppr[:,i], focus_data.ppz[:,i], s=dotsize, color=color)
    plt.axis('equal')
    plt.xlabel('R [m]',fontsize=20)
    plt.ylabel('Z [m]',fontsize=20)
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    # plot iota
    r = np.sqrt(focus_data.ppr[0,:]**2 + focus_data.ppz[0,:]**2)
    plt.figure()
    plt.plot(r, focus_data.iota, marker='s')
    plt.xlabel('R [m]',fontsize=20)
    plt.ylabel(r'$\iota$',fontsize=20)
    plt.xticks(fontsize=16)
    plt.yticks(fontsize=16)
    return

3. Misc.
I will merge this into master after testing. Any comments or suggestions are welcome.

[lazersos]

This seems to work. I just did a test against the FIELDLINES code using a coils file and the two codes seem to qualitatively agree.

[zhucaoxiang]

@lazersos Great, thank you. There are some bugs related the filed-line tracing, in some particular cases. I will push updates and fix them later.