Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Better ways to calculate toroidal flux #44

Open
zhucaoxiang opened this issue Jan 21, 2020 · 18 comments
Open

Better ways to calculate toroidal flux #44

zhucaoxiang opened this issue Jan 21, 2020 · 18 comments
Labels
question

Comments

@zhucaoxiang
Copy link
Collaborator

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

  1. The subroutine of calculating B can be re-used directly. Thus it would be consistent when we changed the bfield.f90 file.
  2. 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.

  1. 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 = 1 to 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.
  2. Can we use another format of current constraint? It is common to specify the overall toroidal current and poloidal current.
  3. 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?
@zhucaoxiang zhucaoxiang mentioned this issue Jan 22, 2020
Closed
@zhucaoxiang
Copy link
Collaborator Author

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

@lazersos
Copy link
Collaborator

@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  ;

@zhucaoxiang
Copy link
Collaborator Author

@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.

@lazersos
Copy link
Collaborator

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

@zhucaoxiang
Copy link
Collaborator Author

@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.

@lazersos
Copy link
Collaborator

@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

@zhucaoxiang
Copy link
Collaborator Author

@lazersos The Bnormal is so high. Are the coils intersecting the plasma?

@lazersos
Copy link
Collaborator

The currents were zero in the TF/VF set. I set them to a finite value and it appears to be running through.

@zhucaoxiang
Copy link
Collaborator Author

When you didn't find errors, you can use xfocus, which is optimized and much faster.

@zhucaoxiang
Copy link
Collaborator Author

@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.

@zhucaoxiang
Copy link
Collaborator Author

@lazersos I just push a commit in symmetry. You should be able to use LM while turning on LC for coil_type=3.

@lazersos
Copy link
Collaborator

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

image

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
Copy link
Collaborator

Of course there's the fact we need to define rho.

@zhucaoxiang
Copy link
Collaborator Author

@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.

@lazersos
Copy link
Collaborator

@zhucaoxiang Would it work to set A_R=Cz/r. Then Bphi=curl(A) = C1/r?

@lazersos
Copy link
Collaborator

@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)

@zhucaoxiang
Copy link
Collaborator Author

@lazersos The answer upvoted looks reasonable and promising.

@zhucaoxiang
Copy link
Collaborator Author

@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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Assignees
No one assigned
Labels
question
Projects
None yet
Milestone
No milestone
Development

No branches or pull requests
2 participants
@lazersos @zhucaoxiang