Vertical Field Coils

[lazersos]

When I try to simulate a Helmholtz pair I find that the current in each coil diverges and the the stellarator coils try to makeup the field. It would be nice if we could include a vertical field as an optimization quantity. Basically one free parameter which could be fixed or varied controlling Bz.

[zhucaoxiang]

@lazersos Rather than changing the source code, there is an easy way to do it. You could include the background field as an target nonzero Bn distribution on the plasma boundary. This could be done either by calculating Bz dot n on the boundary or reading Helmholtz coils, computing Bn, writing updated target Bn by FOCUS. The function of updating the plasma.boundary file is in the Old version. There will be a new file named *.plasma generated, which contains the Bn produced by existing coils. I can also show you when you have time.

[lazersos]

@zhucaoxiang The problem is that I don't know what Bz I want. What I want to do is have the code vary the Bz and the stellarator coil shapes together to target a value of B dot n on the plasma surface. The problem is that the helmholtz pair often drift in their current so that the difference in their currents results in a non-stellarator symmetric change in the coils.

[zhucaoxiang]

@lazersos I think what matters is that the ration of currents in Helmholtz coils and stellarators coils. You can always multiply the currents in both Helmholtz and stellarator coils by a scaling factor, as long as there are no plasma currents presented. How about taking a pair of Helmholtz coils with unit currents (like 1.0 MA), calculating the Bn and using it as the target Bn on the boundary. Afterwards, optimize the shapes and currents in stellarator coils.

[lazersos]

@zhucaoxiang I think this should be expanded to include both a 1/R field and vertical field. Many QAS configurations have net toroidal currents at finite beta so inclusion of vertical field would be useful. Also it's been shown that inclusion of toroidal field coils can help reduce the demands on coil shape for modular coils. Could the 1/R and vertical field be implemented?

[zhucaoxiang]

@lazersos It's on the way. I added one more representation for background magnetic field. The format is something like

 #-----------------           1 ---------------------------
 #coil_type     coil_name
     3    bg_BtBz_01
 # Ic     I    Lc  Bz  (Ic control I; Lc control Bz)
  1  1.000000000000000E+06  0  0.000000000000000E+00

It's already in the branch permanent. But I haven't finished all the implementation. I will keep you informed when it's ready to use. Any suggestions or comments are welcome.

[lazersos]

@zhucaoxiang Found a bug in mapcoil and pushed an update. It looks like everything was implemented in bfield0 but in bfield1 only the vertical field was implemented. Since Bx and By are simply functions of current, wouldn't the derivative just be the same as the bfield0 but without the current term?

[zhucaoxiang]

@lazersos You are correct and that what I have done in the code. If you look at line 147-160 in bnormal.h, I calculate the derivatives of currents by calling bfield0.

Honestly, this branch need deep cleaning and more development to be released. I implemented new coil types, but it's really messy in the code and only Bn is differentiated.

[lazersos]

@zhucaoxiang I tried testing this case but I get errors from Differential Flow, Conjugate Gradient, and Levenberg. Here's the log file:

 ---------------------  FOCUS v0.8.01   ------------------------------
focus   : Begin execution with ncpu =  256
initial : machine_prec   =  2.22045E-16 ; sqrtmachprec   =  1.49012E-08
 -----------INPUT NAMELIST------------------------------------
initial : Read namelist focusin from : QAS.input
        : Read plasma boundary  from : plasma.boundary
        : Read initial coils    from : QAS.focus(Parameters only)
        :         IsQuiet =     -1 ; Output detailed information.
        :     IsSymmetric =      2 ; Plasma boundary and coil periodicity are both enforced.
        :    case_surface =      0 ; Read VMEC-like Fourier harmonics for plasma boundary.
        :      case_coils =      1 ; Using Fourier series as the basic representation.
        :   case_optimize =      1 ; Gradient-based optimizations will be performed.
                          : Differential flow will be used, maxiter=  1000
                          : DF_tausta =  0.00000E+00 ; DF_tauend =  1.00000E+00 ; DF_xtol =  1.00000E-08
        :     IsNormalize =      1 ; Normalization on coil parameters will be performed.
        :    IsNormWeight =      1 ; Normalization on weights for constraints will be performed.
        :    case_bnormal =      1 ; Bn normalized to |B|.
        :     case_length =      2 ; Exponential format of length penalty.
        :   case_postproc =      3 ; Coil evaluations and field-line tracing will be performed.
outputs : HDF5 outputs           are saved in : focus_QAS.h5
outputs : Optimizated coils      are saved in : QAS.focus ; QAS.coils
 -----------Reading surface-----------------------------------
surface : Plasma boundary will be discretized in    128 X    128 elements.
        : Nfou =    137 ; Nfp =      3 ; NBnf =      0 ;
        : Enclosed total plasma volume = 8.14173E+01 m^3 ;
 -----------INITIALIZE COILS----------------------------------
rdcoils : identified      7 unique coils in QAS.focus ;
        :   0 fixed currents ;   0 fixed geometries.
        : Parameter normalizations : Inorm= 3.96973E+06  Gnorm= 1.60372E+00  Mnorm= 1.00000E+00
        : coils will be discretized in    256 segments
focus   : Initialization took  2.44534E-01 S.
 -----------OPTIMIZATIONS-------------------------------------
solvers : Total number of DOF is    278
        : Initial weights are:        bnorm,       bharm,       tflux,       ttlen,       cssep,
        :                       1.00000E+02, 0.00000E+00, 1.00000E+00, 1.00000E+00, 1.00000E+00,
        : target_tflux =  6.55895E+00 ; target_length =  0.00000E+00 ; cssep_factor =  1.00000E+00
 -----------COIL DIAGNOSTICS----------------------------------
diagnos :      Bnormal ; Bmn harmonic ;    tor. flux ;  coil length ;     c-s sep. ; 
        :  2.25046E-02 ;  0.00000E+00 ;  8.39669E-04 ;  1.29974E+08 ;  1.15755E+03 ; 
        : Maximum curvature of all the coils is  :  7.321878671322461E+00 ; at coil   6
        : Average length of the coils is         :  1.614544187228457E+01
        : The minimum coil-coil distance is      :  2.564121606100029E-01
        : The minimum coil-plasma distance is    :  6.298233682289972E-01 ; at coil   1
        : Average relative absolute Bn error is  :  2.005356130777448E-02
        : rescale coil currents with a factor of  1.00004E+00
        : weight_tflux is normalized to  5.12323E+04
        : weight_bnorm is normalized to  4.44331E+03
        : weight_ttlen is normalized to  7.69386E-09
        : weight_cssep is normalized to  8.63891E-04
 ------------- Initial status ------------------------
output  :   iout :         mark ;          chi ;      dE_norm ;      Bnormal ; Bmn harmonic ;    tor. flux ;  coil length ;     c-s sep. ; 
-------------- ERROR ------------------------
length  :      fatal : myid=  1 ; idof .ne. Ndof ; counting error in packing  ;
application called MPI_Abort(MPI_COMM_WORLD, 1) - process 1

I've attached the necessary input files for you to verify the behavior.

test_case.zip

[zhucaoxiang]

@lazersos Thanks for the report. Please pull the branch develop and test it. I just fix this. I tested it with your case and it worked. P.S. It is suggested to use CG. You can check issue #35.

[zhucaoxiang]

@lazersos Please let me know if the commit fixes the problem, such that I can merge it into master.

[lazersos]

@zhucaoxiang It seems to work, isn't crashing the code. So I'd say it's ready for merge into master.