Support reading surface data from an HDF5 file

[zhucaoxiang]

Now users can use an HDF5 file providing surface data.

• To use it, please specifycase_surface = 5 & input_surf = "test.h5". (The option to read Boozer surface is now case_surface=3.)
• In the HDF5 files, you have to provide the following datasets. For 2D arrays, the order must be Fortran-order.
  • Nfp: int, number of toroidal periodicity
  • IsSymmetric: int, stellarator symmetry option (0: no syemmtry; 1: toroidally periodic; 2: periodic and stellarator symmetric)
  • disfac: double precision, the discretization factor for surface integration, usually is 2*pi/Ntheta * 2*pi/Nzeta
  • xsurf: double precision (Ntheta*Nzeta), x coordinates of the surface elements
  • ysurf: double precision (Ntheta*Nzeta), y coordinates of the surface elements
  • zsurf: double precision (Ntheta*Nzeta), z coordinates of the surface elements
  • nx: double precision (Ntheta*Nzeta), x coordinates of the surface unit normal
  • ny: double precision (Ntheta*Nzeta), y coordinates of the surface unit normal
  • nz: double precision (Ntheta*Nzeta), z coordinates of the surface unit normal
  • nn: double precision (Ntheta*Nzeta), the surface jacobian (normal vector magnitude)
  • plas_Bn: double precision (Ntheta*Nzeta), target Bn (FOCUS uses Bn_coil - Bn_target
• You can also use the FOCUS-saved HDF5 file to provide surface information.

There are also some limitations.

• The dr/dtheta and dr/dzeta components are not read. So torflux is not computed correctly. One can fix or constrain coil currents, instead.
• When initializing circular coils or providing an initial guess for the magnetic axis, the subroutinesurfcood is called. It reads (theta, zeta) and returns (R, Z). When reading the HDF files, Fourier coefficients for the surface are unknown. So surfcoord now returns the nearest point in the given surface data. In short, initial circular coils and initial magnetic axis will be slightly different.

[zhucaoxiang]

Here is a verification case. test_fourier uses case_surface = 0 and performs a simple optimization from circular coils. test uses case_surface = 5, reads focus_test_fourier.h5, and performs the same optimization.

Here is part of the results for test_fourier.

output  :    100 :  3.26495E+01 ;  3.58817E-03 ;  1.17884E-03 ;  1.14991E-03 ;  6.21720E-04 ;
output  :    101 :  3.30023E+01 ;  3.58787E-03 ;  1.25395E-03 ;  1.14995E-03 ;  6.21148E-04 ;
congrad : status=2: total iterations exceeded maxit.
congrad : Computation using conjugate gradient finished.
        : CG takes   3.300243400000000E+01 seconds;
focus   : Optimization took  3.33170E+01 S.
 -----------POST-PROCESSING-----------------------------------
 -----------COIL DIAGNOSTICS----------------------------------
diagnos : Individual cost function values:
        :              Bnormal =  1.14994787E-03
        :       B_mn harmonics =  0.00000000E+00
        :        Toroidal flux =  0.00000000E+00
        : Total current penalt =  0.00000000E+00
        :          Coil length =  6.21148232E-04
        :       Coil curvature =  0.00000000E+00
        :         Coil torsion =  0.00000000E+00
        :    Nissin complexity =  0.00000000E+00
        : Coil-coil separation =  0.00000000E+00
        : Coil-surf separation =  0.00000000E+00
        :   Stochastic Bnormal =  0.00000000E+00
        :         Straight-out =  0.00000000E+00
        : Average length of the coils is         :  1.033895324843059E+01
        : Maximum curvature of all the coils is  :  5.562765354252051E+00 ; at coil   2
        : Average curvature of the coils is      :  1.119370574850224E+00
        : Average torsion of the coils is        :  7.798044489677238E-01
        : The minimum coil-coil distance is      :  2.342042451264976E-01 ; at coils  2 ,  1
        : The minimum coil-plasma distance is    :  4.615120032344702E-01 ; at coil   1
        : Ave. relative absolute Bn error |Bn|/B :  1.07866E-02; max(|Bn|)= 5.35223E-02
        : Surface area normalized Bn error int(|Bn|/B*ds)/A :   1.007061748653700E-02

Here is part of the results for test. They are identical (Computational time is slightly different, due to machine performance fluctuations).

output  :    100 :  3.21527E+01 ;  3.58817E-03 ;  1.17884E-03 ;  1.14991E-03 ;  6.21720E-04 ;
output  :    101 :  3.24544E+01 ;  3.58787E-03 ;  1.25395E-03 ;  1.14995E-03 ;  6.21148E-04 ;
congrad : status=2: total iterations exceeded maxit.
congrad : Computation using conjugate gradient finished.
        : CG takes   3.245481700000001E+01 seconds;
focus   : Optimization took  3.27795E+01 S.
 -----------POST-PROCESSING-----------------------------------
 -----------COIL DIAGNOSTICS----------------------------------
diagnos : Individual cost function values:
        :              Bnormal =  1.14994787E-03
        :       B_mn harmonics =  0.00000000E+00
        :        Toroidal flux =  0.00000000E+00
        : Total current penalt =  0.00000000E+00
        :          Coil length =  6.21148232E-04
        :       Coil curvature =  0.00000000E+00
        :         Coil torsion =  0.00000000E+00
        :    Nissin complexity =  0.00000000E+00
        : Coil-coil separation =  0.00000000E+00
        : Coil-surf separation =  0.00000000E+00
        :   Stochastic Bnormal =  0.00000000E+00
        :         Straight-out =  0.00000000E+00
        : Average length of the coils is         :  1.033895324835169E+01
        : Maximum curvature of all the coils is  :  5.562765369656166E+00 ; at coil   2
        : Average curvature of the coils is      :  1.119370574856432E+00
        : Average torsion of the coils is        :  7.798044486602599E-01
        : The minimum coil-coil distance is      :  2.342042450923805E-01 ; at coils  2 ,  1
        : The minimum coil-plasma distance is    :  4.615120032998292E-01 ; at coil   1
        : Ave. relative absolute Bn error |Bn|/B :  1.07866E-02; max(|Bn|)= 5.35223E-02
        : Surface area normalized Bn error int(|Bn|/B*ds)/A :   1.007061748492264E-02

hdf5test.zip