PrincetonUniversity
Better ways to calculate toroidal flux
There is a cost function on the toroidal flux in FOCUS. Instead of calculating B \cdot ds, FOCUS uses the Stokes theorem and calculates A \cdot dl. Details can be found in torflux.pdf.
This works fine with normal coils. Recently, I found that it might be better to use B \cdot ds, because
- The subroutine of calculating B can be re-used directly. Thus it would be consistent when we changed the bfield.f90 file.
- The vector potential of an infinitely long current wire is actually infinite Eq. 9.3.5, which means using A-field doesn't work for the toroidal field generated by a current in the center of the torus. (Maybe I am wrong?)
To use B \cdot ds, one might select the phi=0 cut. I realize that it requires some non-trivial effort to get the toroidal cut. Basically, one needs to check if an arbitrary point on the XZ plane is inside or outside the plasma. This could be solved by counting the winding number or cross number wikipedia.
Here are some questions that we can discuss.
- Do we need to use the toroidal flux constraint? It was imported to avoid trivial solutions of all currents going to zero. By fixing the coil current, specifying non-zero target Bn (e.g. from plasma), or using
case_bnormal = 1to optimized B \cdot n / |B|, one can also avoid such a trivial solution. The other good side of using the toroidal flux is to incorporate with equilibrium codes, like VMEC and SPEC. In such codes, the edge toroidal flux is always specified. - Can we use another format of current constraint? It is common to specify the overall toroidal current and poloidal current.
- If we are going to use the toroidal flux and use the B \cdot ds, is there an easier way to integrate B \cdot ds?
Right now, torflux will return an error when using coil_type other than 1, because of the limitation on using vector potential. If you want to use other coil_type, do NOT use toroidal flux constraint. Instead, you can use case_Bnormal = 1 or just simply fix one (or all) current.
@lazersos
@zhucaoxiang I attempted to comment out the two lines you suggested and now I get this error:
-----------INITIALIZE COILS----------------------------------
rdcoils : identified 2 unique coils in QASDEX.focus ;
: 0 fixed currents ; 0 fixed geometries.
: Parameter normalizations : Inorm= 7.07107E+05 Gnorm= 4.58291E+00 Mnorm= 1.00000E+00
: coils will be discretized in 256 segments
focus : Initialization took 1.43436E+00 S.
-----------OPTIMIZATIONS-------------------------------------
solvers : Total number of DOF is 24
: Initial weights are: bnorm, bharm, tflux, ttlen, cssep,
: 1.00000E+02, 0.00000E+00, 0.00000E+00, 1.00000E+00, 1.00000E+00,
: target_tflux = 0.00000E+00 ; target_length = 0.00000E+00 ; cssep_factor = 4.00000E+00
datalloc: total number of cost functions for L-M is 8196
-----------COIL DIAGNOSTICS----------------------------------
costfun : Reset target toroidal flux to 1.114269080647567E+00
-------------- ERROR ------------------------
length : fatal : myid= 0 ; ivec == mttlen ; Errors in counting ivec for L-M ;
@lazersos Good catch. I haven't used LM for a long while and didn't test it before my commit. Could you please use CG temporarily? I will check what is going on in LM.
My personal feeling for LM in FOCUS is that it is more violent than CG. Maybe more efficient but also tends to provide over-optimized results.
New error with the CG
-----------INITIALIZE COILS----------------------------------
rdcoils : identified 2 unique coils in QASDEX.focus ;
: 0 fixed currents ; 0 fixed geometries.
: Parameter normalizations : Inorm= 7.07107E+05 Gnorm= 4.58291E+00 Mnorm= 1.00000E+00
: coils will be discretized in 256 segments
focus : Initialization took 1.30164E+00 S.
-----------OPTIMIZATIONS-------------------------------------
solvers : Total number of DOF is 24
: Initial weights are: bnorm, bharm, tflux, ttlen, cssep,
: 1.00000E+02, 0.00000E+00, 0.00000E+00, 1.00000E+00, 1.00000E+00,
: target_tflux = 0.00000E+00 ; target_length = 0.00000E+00 ; cssep_factor = 4.00000E+00
-----------COIL DIAGNOSTICS----------------------------------
costfun : Reset target toroidal flux to 1.114269080647567E+00
diagnos : Bnormal ; Bmn harmonic ; tor. flux ; coil length ; c-s sep. ;
: 2.71096E+03 ; 0.00000E+00 ; 6.95376E-02 ; 4.26837E+40 ; 2.93314E+00 ;
: Maximum curvature of 1-th coil is : 6.666666666666668E-02
: Maximum curvature of all the coils is : 6.666666666666668E-02 ; at coil 1
: Average length of the coils is : 4.712389111518860E+01
: distance between 001-th and 002-th coil (ip=01, is=0) is : 1.000000000000000E+03
: distance between 001-th and 002-th coil (ip=02, is=0) is : 1.000000000000000E+03
: The minimum coil-coil distance is : 1.000000000000000E+03
: The minimum coil-plasma distance is : 8.068390097051161E+00 ; at coil 1
: Ave. relative absolute Bn error |Bn|/B : 3.15238E+00; max(|Bn|)= 2.57851E-01
: Surface area normalized Bn error int(|Bn|/B*ds)/A : 2.787976862161040E+00
: weight_bnorm is normalized to 3.68873E-02
: weight_ttlen is normalized to 2.34282E-41
warning : weight_ttlen < machine_precision, ttlen will not be used.
: weight_cssep is normalized to 3.40931E-01
------------- Initial status ------------------------
output : iout : mark ; chi ; dE_norm ; Bnormal ; Bmn harmonic ; tor. flux ; coil length ; c-s sep. ;
output : 1 : 0.00000E+00 ; 1.01000E+02 ; NaN ; 2.71096E+03 ; 0.00000E+00 ; 6.95376E-02 ; 4.26837E+40 ; 2.93314E+00 ;
forrtl: severe (174): SIGSEGV, segmentation fault occurred
Image PC Routine Line Source
dfocus 00000000005690CD Unknown Unknown Unknown
libpthread-2.22.s 00002AFAFB52EC70 Unknown Unknown Unknown
dfocus 00000000004F1FFC saving_ 4751 saving_m.F90
dfocus 00000000004BEF75 output_ 585 solvers_m.F90
dfocus 00000000004B7353 solvers_ 90 solvers_m.F90
dfocus 000000000051C317 MAIN__ 102 focus_m.F90
dfocus 000000000040B71E Unknown Unknown Unknown
libc-2.22.so 00002AFAFD0B3765 __libc_start_main Unknown Unknown
dfocus 000000000040B629 Unknown Unknown Unknown
@lazersos It was caused by writing MAKEGIRD format coils for coil_type=3 that doesn't have x,y,z data. I have fixed it. It will skip coil_type .ne. 1.
@zhucaoxiang I think there are still issues
-----------COIL DIAGNOSTICS----------------------------------
costfun : Reset target toroidal flux to 1.114269080647567E+00
diagnos : Bnormal ; Bmn harmonic ; tor. flux ; coil length ; c-s sep. ;
: 2.71096E+03 ; 0.00000E+00 ; 6.95376E-02 ; 4.26837E+40 ; 2.93314E+00 ;
: Maximum curvature of 1-th coil is : 6.666666666666668E-02
: Maximum curvature of all the coils is : 6.666666666666668E-02 ; at coil 1
: Average length of the coils is : 4.712389111518860E+01
: distance between 001-th and 002-th coil (ip=01, is=0) is : 1.000000000000000E+03
: distance between 001-th and 002-th coil (ip=02, is=0) is : 1.000000000000000E+03
: The minimum coil-coil distance is : 1.000000000000000E+03
: The minimum coil-plasma distance is : 8.068390097051161E+00 ; at coil 1
: Ave. relative absolute Bn error |Bn|/B : 3.15238E+00; max(|Bn|)= 2.57851E-01
: Surface area normalized Bn error int(|Bn|/B*ds)/A : 2.787976862161040E+00
: weight_bnorm is normalized to 3.68873E-02
: weight_ttlen is normalized to 2.34282E-41
warning : weight_ttlen < machine_precision, ttlen will not be used.
: weight_cssep is normalized to 3.40931E-01
------------- Initial status ------------------------
output : iout : mark ; chi ; dE_norm ; Bnormal ; Bmn harmonic ; tor. flux ; coil length ; c-s sep. ;
output : 1 : 0.00000E+00 ; 1.01000E+02 ; NaN ; 2.71096E+03 ; 0.00000E+00 ; 6.95376E-02 ; 4.26837E+40 ; 2.93314E+00 ;
------------- Nonlinear Conjugate Gradient (CG) -----
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
iter: 0 f= 0.101000E+03 gnorm= 0.871790E+03 AWolfe= F
QuadOK: F initial a: 0.750875E-04 f0: 0.101000E+03 dphi NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
secant, a: 0.000000E+00 b: 0.750875E-04 da: NaN db: NaN
@lazersos The Bnormal is so high. Are the coils intersecting the plasma?
The currents were zero in the TF/VF set. I set them to a finite value and it appears to be running through.
When you didn't find errors, you can use xfocus, which is optimized and much faster.
@lazersos coil_type=3 was added later and there are some incompatibilities that I didn't catch when I merged, like LM. Here is my suggestion.
For the background field, try to use Bz = 1.0 and Lc=0 to fix the vertical field. Set I = 1E6 (or a finite number and Ic = 1 to vary the central current. For other normal coils (type=1) you can turn both Ic and Lc on. For your cost functions, you don't have to use toroidal flux, because there is a 1T vertical field to cancel, coil current cannot be zero.
@lazersos I just push a commit in symmetry. You should be able to use LM while turning on LC for coil_type=3.
I had a thought of on this problem. You're right that the vector potential for a current carrying wire is infinite but for a solenoid it's well defined. So if we use
We just need to change nI to I since nI really just represents the total current. It's approximate but since in the end we'll need to design actual TF coils this is probably just as reasonable as assuming an infinite wire on axis.
@lazersos I was thinking to use B dot n ds to calculate the toroidal flux. We can use the phi=0 cross-section where the normal vector is the y-axis. The only headache is that we have to determine which surface element on the x-z plane is inside the boundary. It is workable, but I didn;t find an elegant way to do this.
@zhucaoxiang or you could use the logarithmic form (https://physics.stackexchange.com/questions/222343/magnetic-vector-potential-of-an-infinite-wire). A_Z=C*ln(r)
@lazersos I have implemented the formula in toroidal flux. But I didn't fully verify and validate it. Please let me know if you found any bug when using it.