search

PrincetonUniversity / SPEC

The Stepped-Pressure Equilibrium Code, an advanced MRxMHD equilibrium solver.
https://princetonuniversity.github.io/SPEC/
GNU General Public License v3.0
24 stars 4 forks source link

VMEC-SPEC benchmark #70

Closed zhucaoxiang closed 4 years ago

zhucaoxiang commented 5 years ago

@zhisong @jloizu CC: @SRHudson I am working on SPEC-VMEC benchmark. Right now, to get free-boundary SPEC converged, you need to provide the exact Vns from coils (calculated by Biot-Savart law) and a good initial guess for Bns from plasma current on the computational boundary.

Zhisong did the following to get it converged.

But there is a more friendly way, which is SPEC calculates the initial guess for Bns from fixed-boundary run. This will normally be a good guess. I added the following line to let SPEC run fixed-boundary, do virtual casing and update Bns at the first step of free-boundary calculation. https://github.com/PrincetonUniversity/SPEC/blob/5e0e4c775561ffa5a558317a7d2fbd3c2e2b8026/xspech.h#L172-L176

But somehow, after the fixed boundary run, the virtual casing is not converged, as shown below.

o smal;        
casing :     179.39 : myid= 14 ; [x,y,z]=[  5.39E+00 ,  0.00E+00 ,  1.14E+00 ]; gBn=  4.2676E-01 , err=  3.E-03 ; ifail=  1 ; min/max calls=    16777215    33554432 ; maxpts too smal;        
casing :     179.52 : myid=  5 ; [x,y,z]=[  5.75E+00 ,  0.00E+00 ,  4.37E-01 ]; gBn= -1.3139E-01 , err=  2.E-02 ; ifail=  1 ; min/max calls=    16777215    33554432 ; maxpts too smal;        
casing :     179.69 : myid=  0 ; [x,y,z]=[  5.80E+00 ,  0.00E+00 ,  0.00E+00 ]; gBn= -7.3118E-10 , err=  3.E-03 ; ifail=  1 ; min/max calls=    16777215    33554432 ; maxpts too smal;        
casing :     180.36 : myid=  8 ; [x,y,z]=[  5.66E+00 ,  0.00E+00 ,  6.89E-01 ]; gBn=  1.5996E-01 , err=  3.E-03 ; ifail=  1 ; min/max calls=    16777215    33554432 ; maxpts too smal;        
casing :     180.51 : myid=  2 ; [x,y,z]=[  5.79E+00 ,  0.00E+00 ,  1.76E-01 ]; gBn=  3.5361E-02 , err=  3.E-03 ; ifail=  1 ; min/max calls=    16777215    33554432 ; maxpts too small;

I need your help, as you are more familiar with SPEC source code. Please pull the branch free_bound, compile it and run the example in InputFiles/Verification/forcefree/solovev/solovev_fb_1r.sp. In this input file, I provide the required Vns from FOCUS and zero for Bns, which is the general case for users. I set gBnbld=0 to using virtual casing from fixed-boundary run, but probably we need to reset it to a larger number after Bns was obtained.

zhisong commented 5 years ago

I believe in casing.h, the virtual casing integral takes the field on the the inner interface of the vacuum region. Your modification to run SPEC without vacuum region will not work, since it does not provide solution to the vacuum region.

The reason why I ran virtual casing 30 times is that, the field on the inner interface of the vacuum region and Bns depend on each other, and thus one needs an iteration process(Picard) to get both at the same time. I don't see a way around at the moment.

zhucaoxiang commented 5 years ago

@zhisong

the inner interface of the vacuum region

Are you talking about the plasma boundary or the computational boundary? In principle, this method should work. Once you finish a fixed-boundary run, you can use virtual casing to calculate the magnetic field outside the plasma, here we use the computation boundary. I don't think the solution in the vacuum region is relevant.

Before calling Bnormal for virtual casing, I recover Mvol=Nvol+Lfreebound to make sure it is using the computational boundary.

There is actually one thing I forgot. @SRHudson mentioned there is a flag in virtual casing calculation to specify using the interior or exterior surface current. There might be difference if existing a sheet current on the plasma boundary. But I didn't see anything

zhisong commented 5 years ago

Let me say it in another way. Currently virtual casing is using the field on the exterior of the plasma boundary. Therefore, you need to have the vacuum solution ready. Otherwise, you don't know what the field is.

In general, you can use the field on the interior of the plasma boundary, if you don't care about the surface current on the plasma boundary. However, this option is not ready-to-use. One will need to modify the source code. I did hack the source to use the interior when we performed the wire benchmark. But then Joaquim changed it back.

If, however, one does care about that surface current, he/she will need to run virtual casing 30 times as I did, although one can use the interior field to get a better initial guess and reduce the number of iteration to, say, 10.

zhucaoxiang commented 5 years ago

@zhisong That is what I mean. For the first attempt, we can just use the interior to provide a considerably good initial guess. It might not be accurate, but it will be close enough as long as there is no large sheet current on the plasma boundary. I don't see any reason that SPEC needs 10 iterations. One fixed-boundary run should be enough.

Here is my idea. For a general free-boundary case, the user provides Vns from external coils, plasma boundary and some profiles (tflux, iota, etc). If Bns are zero (not provided), for the first iteration, SPEC will

  1. run in fixed-boundary mode with the given plasma boundary;
  2. call virtual casing to calculate Bns on the computational boundary using the interior surface current;
  3. update Bns with gBnbld=0;
  4. begin regular free-boundary iterations with updated Bns.
zhucaoxiang commented 5 years ago

@zhisong I added the following lines in casing.h trying to use the interior. https://github.com/PrincetonUniversity/SPEC/blob/9bffb8b1a7e412a762e8785bb18b704fcc9d7746/casing.h#L396-L406

But there are two more remaining issues.

  1. Virtual casing precision for nfreebounditerations>1 cannot be as low as 1E-10. It is something like 1E-5.
  2. abs(det).lt.small error. 
    fcn2   :      73.78 :        29  1 ; |f|= 4.75493E-15 ; time=      1.07s ; log|BB|e=-15.63-15.96-17.12-17.52-17.21-17.16
fcn2   :            :              ;                                     ; log|II|o=-14.39-14.48-14.44-14.19-14.17-14.42
newton :            :
newton :      73.78 : finished ; success        ; ic05p*f= 1 ; its=     29 ,   1 ;
xspech :      74.81 : #freeits=  0 ; |f|= 4.75493E-15 ; time=      1.00s ; log|BB|e=-15.63-15.96-17.12-17.52-17.21-17.16
xspech :            :              ;                                     ; log|II|o=-14.39-14.48-14.44-14.19-14.17-14.42
xspech :            : 
xspech :     104.33 : nfreeboundaryiterations =      0 /  00020 ; gBntol = 1.0E-10 ; bnserr = 1.05012E-03 ; bnorml time =     29.52s ;
xspech :            : 
nfreeboundaryiterations=           1
ma02aa :     106.86 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  4 hel= -3.2379E-03 mu=  0.0000E+00 dpflux= -1.1886E-02 time=      1.5 ; bad progress     ; F=  2.E-01 -1.E-03
newton :     107.03 :         0  0 ; |f|= 9.34332E-03 ; time=      2.68s ; log|BB|e=-15.63-15.96-17.12-17.52-17.21-17.15 -2.42
newton :            :              ;                                     ; log|II|o=-14.39-14.48-14.44-14.19-14.17-14.42 -1.78
ma02aa :     109.25 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -1.2100E-01 mu=  0.0000E+00 dpflux= -4.8388E-02 time=      1.1 ; bad progress     ; F=  4.E-02 -2.E-06
ma02aa :     110.98 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -1.2048E-01 mu=  0.0000E+00 dpflux= -4.8310E-02 time=      0.6 ; bad progress     ; F=  4.E-02 -5.E-05
fcn2   :     150.39 :         1  1 ; |f|= 9.69386E-03 ; time=     46.05s ; log|BB|e=-15.63-15.96-17.12-17.52-17.21-17.15 -2.06
fcn2   :            :              ;                                     ; log|II|o=-14.39-14.48-14.44-14.19-14.17-14.42 -1.78
ma02aa :     152.63 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -9.3385E-02 mu=  0.0000E+00 dpflux= -3.9747E-02 time=      1.3 ; bad progress     ; F=  5.E-02 -1.E-03
fcn2   :     152.88 :         2  1 ; |f|= 8.31782E-03 ; time=      2.48s ; log|BB|e=-15.85-16.33-16.78-17.21 -5.66 -5.40 -2.27
fcn2   :            :              ;                                     ; log|II|o=-14.53-14.63-14.60-14.34 -3.31 -2.79 -1.83
fcn2   :     154.32 :         3  1 ; |f|= 6.20851E-03 ; time=      1.44s ; log|BB|e=-16.51-16.22 -6.35 -5.80 -4.77 -4.45 -2.83
fcn2   :            :              ;                                     ; log|II|o=-15.36-15.41 -4.48 -3.44 -3.39 -2.45 -1.97
ma02aa :     157.27 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -2.5396E-01 mu=  0.0000E+00 dpflux= -6.2672E-02 time=      1.9 ; bad progress     ; F= -7.E-02  8.E-03
fcn2   :     157.56 :         4  1 ; |f|= 4.89864E-03 ; time=      3.24s ; log|BB|e= -6.93 -6.21 -5.33 -4.63 -5.09 -4.20 -2.49
fcn2   :            :              ;                                     ; log|II|o= -5.94 -4.40 -3.72 -3.23 -2.66 -2.73 -2.09
ma02aa :     160.49 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -1.4774E-02 mu=  0.0000E+00 dpflux=  7.0837E-02 time=      2.1 ; bad progress     ; F=  1.E-01 -2.E-04
fcn2   :     160.73 :         5  1 ; |f|= 6.38632E-03 ; time=      3.17s ; log|BB|e= -5.86 -5.54 -5.82 -4.72 -5.18 -4.17 -1.82
fcn2   :            :              ;                                     ; log|II|o= -5.37 -4.87 -3.58 -3.52 -3.36 -2.69 -2.10
ma02aa :     162.24 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -1.6545E+05 mu=  0.0000E+00 dpflux= -3.2104E+04 time=      0.7 ; bad progress     ; F=  2.E-01  8.E-03
ma02aa :     164.29 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel= -1.2369E+13 mu=  0.0000E+00 dpflux= -6.6373E+12 time=      0.8 ; bad progress     ; F=  2.E-01 -7.E-03
dforce :      fatal : myid=  7 ; abs(det).lt.small ; error computing derivatives of dtflux & dpflux wrt geometry at fixed Itor and Gpol ;
zhisong commented 5 years ago

I understand your point. In principle it should work.

However, I am unsure about whether simply set Mvol=Nvol is a clean way to do so. One needs to go to the code lines to check if the logic involves some other variables such as Lfreebound.

Alternatively, if we don't care about user-friendliness, we can run a sequence of three input files: first with Lfreebound=false and Lfindzero=2, then Lfreebound=true, Lvirtualcasing=inside and Lfindzero=0, finally with Lfreebound=true, Lvirtualcasing=outside and Lfindzero=2.

jloizu commented 5 years ago

Hello,

I am catching up with the discussion. I agree that in order to get a quick (without any iteration) guess for the Bns on the computational boundary we could do one virtual casing evaluation taking the inner field and using gBnbld=0. But we do could do this directly in free-boundary mode, i.e. no need of doing a fixed-boundary run first.

What I mean is:

  1. run with Lfindzero=0, Lvirtualcasing=inside, gBnbld=0, mfreeits=1 (using Bns=0 as guess)
  2. run with Lfindzero=2, Lvirtualcasing=outside, gBnbld<1, mfreeits>1 (using Bns from previous run)

which is what Zhisong proposed in the last comment, but omitting the first step, since the free-boundary run at fixed geometry (my step 1) already is doing a fixed-boundary calculation inside the plasma boundary, and immediately after the virtual casing would called with the appropriate value for the field inside the plasma boundary).

Does this makes sense? If so, It would be easy to make this change and then one would need to run twice the code for each equilibrium. Then we can make it such that it does this internally.

zhisong commented 5 years ago

@jloizu your first step will work if the internal interfaces are initialized very well, say extracted from VMEC. Otherwise we will need do a fixed boundary run and use the converaged result to perform virtual casing.

jloizu commented 5 years ago

Ok I agree. Sounds ugly thought to need to run 3 input files...shall we try first like that and then if it works we make SPEC do this internally?

zhisong commented 5 years ago

I agree with @jloizu . We first make sure that the procedure works, then deal with user-friendliness. @zhucaoxiang does the proposed sequence-of-three-file method work on on your end?

zhucaoxiang commented 5 years ago

@zhisong @jloizu Seems this is the right approach. I will test it.

zhucaoxiang commented 5 years ago

@zhisong @jloizu I assume Lvirtualcasing=inside is what I did as the following https://github.com/PrincetonUniversity/SPEC/blob/9bffb8b1a7e412a762e8785bb18b704fcc9d7746/casing.h#L396-L400

Somehow, I did Lfindzero=0, Lvirtualcasing=inside, gBnbld=0, mfreeits=1 run with first_free_bound=True and virtual casing is still not converging.

casing :      43.23 : myid= 10 ; [x,y,z]=[  5.59E+00 ,  0.00E+00 ,  8.49E-01 ]; gBn= -8.2066E-02 , err=  1.E-03 ; ifail=  1 ; min/max calls=           8    16777216 ; maxpts too smal;        
casing :      43.34 : myid=  5 ; [x,y,z]=[  5.75E+00 ,  0.00E+00 ,  4.37E-01 ]; gBn=  4.1831E-02 , err=  4.E-03 ; ifail=  1 ; min/max calls=           8    16777216 ; maxpts too smal;        
casing :      43.39 : myid=  7 ; [x,y,z]=[  5.69E+00 ,  0.00E+00 ,  6.06E-01 ]; gBn= -1.3774E-01 , err=  1.E-03 ; ifail=  1 ; min/max calls=           8    16777216 ; maxpts too smal;        
casing :      43.48 : myid=  8 ; [x,y,z]=[  5.66E+00 ,  0.00E+00 ,  6.89E-01 ]; gBn= -5.2811E-02 , err=  1.E-03 ; ifail=  1 ; min/max calls=           8    16777216 ; maxpts too smal;

Anyone has some hints?

jloizu commented 5 years ago

Isn't just that you have put some value for the tolerance gBntol = 1.000000000000000E-10 ? This measures the difference between the obtained (unblended) Bns and the previous one (which in your case is Bns=0). So the virtual casing is not converged in that sense but we don't care. Does this make sense? You can try setting gBntol to a large number if the non-convergence does not allow finishing the job.

zhucaoxiang commented 5 years ago

@jloizu This is determined by vcasingtol and I don't think it is related to the Picard iteration. To my understanding this is something like the tolerance for the surface integral. I definitely can ignore the tolerance, although I didn't quite understand why it doesn't converge.

zhucaoxiang commented 5 years ago

Once I relaxed vcasingtol, I got Bns updated and for the normal free-boundary run with gBnbld=0.8. I got the following error.

preset :       0.00 : LBsequad= F , LBnewton= F , LBlinear= T ;
preset :            : 
preset :       0.01 : Nquad=  -1 ; mn=   17 ; NGdof=   231 ; NAdof=   341,   358,   358,   358,   358,   358,   358,   426,
preset :            : 
preset :       0.01 : Nt=   128 ; Nz=     1 ; Ntz=      128 ;
ma02aa :       1.86 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel=  1.8681E+06 mu=  0.0000E+00 dpflux=  5.7913E+05 time=      1.0 ; bad progress     ; F=  2.E-01  9.E-07
newton :       1.87 :         0  0 ; |f|= 5.59014E+11 ; time=      1.83s ; log|BB|e= -5.12 -4.97 -4.83 -4.68 -4.53 -4.38 12.27
newton :            :              ;                                     ; log|II|o= -4.27 -3.45 -3.03 -2.82 -2.50 -2.25 -1.81
ma02aa :       3.16 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel=  3.3592E+12 mu=  0.0000E+00 dpflux=  5.7913E+05 time=      0.5 ; bad progress     ; F=  2.E-01  9.E-07
ma02aa :       4.50 : myid=  7 ; lvol=  8 ; Linear : ihybrj =  5 hel=  3.3592E+12 mu=  0.0000E+00 dpflux=  5.7913E+05 time=      0.5 ; bad progress     ; F=  2.E-01  9.E-07
dforce :      fatal : myid=  7 ; abs(det).lt.small ; error computing derivatives of dtflux & dpflux wrt geometry at fixed Itor and Gpol ;

My guess is that the calculated Bns is not even close, as the virtual casing is finished with relatively low precision. I will try to enforce vcasingtol=1E-6 and see it can be reached.

jloizu commented 5 years ago

Just as a reference, for the cases we constructed to verify the virtual casing algorithm (the so-called wire model), we had

 vcasingeps  =   1.000000000000000E-12
 vcasingtol  =   1.000000000000000E-08
 vcasingits  =         8
 vcasingper = 1
zhucaoxiang commented 5 years ago

@jloizu Thanks. I just set vcasingtol = 1.0E-6 and I am not optimistic. It is already 134 millions calls. I guess there is some places we missed.

casing :    1222.77 : myid=  6 ; [x,y,z]=[  5.72E+00 ,  0.00E+00 ,  5.23E-01 ]; gBn=  3.7497E-01 , err=  5.E-05 ; ifail=  1 ; min/max calls=   134217655   268435456 ; maxpts too smal;        
casing :    1227.54 : myid=  7 ; [x,y,z]=[  5.69E+00 ,  0.00E+00 ,  6.06E-01 ]; gBn= -8.0345E-01 , err=  5.E-05 ; ifail=  1 ; min/max calls=   134217655   268435456 ; maxpts too smal;        
casing :    1227.68 : myid=  1 ; [x,y,z]=[  5.80E+00 ,  0.00E+00 ,  8.83E-02 ]; gBn=  1.3596E+00 , err=  5.E-05 ; ifail=  1 ; min/max calls=   134217655   268435456 ; maxpts too smal;        
casing :    1230.45 : myid=  5 ; [x,y,z]=[  5.75E+00 ,  0.00E+00 ,  4.37E-01 ]; gBn=  1.1416E+00 , err=  1.E-03 ; ifail=  1 ; min/max calls=   134217655   268435456 ; maxpts too smal;        
casing :    1235.04 : myid= 11 ; [x,y,z]=[  5.54E+00 ,  0.00E+00 ,  9.25E-01 ]; gBn= -1.0110E-01 , err=  8.E-04 ; ifail=  1 ; min/max calls=   134217655   268435456 ; maxpts too smal;        
casing :    1237.93 : myid=  0 ; [x,y,z]=[  5.80E+00 ,  0.00E+00 ,  0.00E+00 ]; gBn= -8.5542E-15 , err=  5.E-05 ; ifail=  1 ; min/max calls=   134217655   268435456 ; maxpts too smal;
zhucaoxiang commented 5 years ago

Here is the input file I am using.

&physicslist
 Igeometry   =         3
 Istellsym   =         1
 Lfreebound  =         1
 phiedge     =   1.000000000000000E+00
 curtor      =  -1.190644396046887E-01
 curpol      =   4.881440045660601E+00
 gamma       =   0.000000000000000E+00
 Nfp         =         1
 Nvol        =         7
 Mpol        =        16
 Ntor        =         0
 Lrad        =                       8                      8                      8                      8                      8                      8                      8                     10
 tflux       =   2.040816326530612E-02  8.163265306122448E-02  1.836734693877551E-01  3.265306122448979E-01  5.102040816326531E-01  7.346938775510207E-01  1.000000000000000E+00  2.091907318303254E+00
 pflux       =   0.000000000000000E+00 -2.884269637968015E-02 -7.440688655913977E-02 -1.329405852326906E-01 -1.991894828261694E-01 -2.664050791017344E-01 -3.264076564461609E-01 -5.032315890880776E-01
 helicity    =  -2.016178262074368E-04  1.559429589793997E-03  1.559429589793997E-03  1.559429589793997E-03  3.042398287639578E-04  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
 pscale      =   1.554767792848010E-07
 Ladiabatic  =         0
 pressure    =   1.000000000000000E+00  9.587628865979375E-01  8.762886597938135E-01  7.525773195876289E-01  5.876288659793815E-01  3.814432989690719E-01  1.340206185567007E-01  0.000000000000000E+00
 adiabatic   =   1.000000000000000E+00  9.587628865979375E-01  8.762886597938135E-01  7.525773195876289E-01  5.876288659793815E-01  3.814432989690719E-01  1.340206185567007E-01  0.000000000000000E+00
 mu          =  -2.500197260488832E-01 -2.392436660969981E-01 -2.175753295323239E-01 -1.847587151649739E-01 -1.404239839402925E-01 -8.458267058227775E-02 -2.077978848646100E-02  0.000000000000000E+00
 Lconstraint =         1
 pl          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 ql          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 pr          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 qr          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 iota        =   0.000000000000000E+00 -4.838775510204081E-01 -4.655102040816324E-01 -4.348979591836737E-01 -3.920408163265306E-01 -3.369387755102044E-01 -2.695918367346938E-01 -1.900000000000000E-01  0.000000000000000E+00
 lp          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 lq          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 rp          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 rq          =                       0                      0                      0                      0                      0                      0                      0                      0                      0
 oita        =   0.000000000000000E+00 -4.838775510204081E-01 -4.655102040816324E-01 -4.348979591836737E-01 -3.920408163265306E-01 -3.369387755102044E-01 -2.695918367346938E-01 -1.900000000000000E-01  0.000000000000000E+00
 mupftol     =   1.000000000000000E-12
 mupfits     =       128
 Rac         =   4.095909789607840E+00
 Zas         =   0.000000000000000E+00
 Ras         =   0.000000000000000E+00
 Zac         =   0.000000000000000E+00
Rbc(0,0)    =  4.040456517412868E+00 Zbs(0,0)    =  0.000000000000000E+00 Rbs(0,0)    =  0.000000000000000E+00 Zbc(0,0)    =  0.000000000000000E+00
Rbc(0,1)    =  1.037511255072165E+00 Zbs(0,1)    = -1.562461763819107E+00 Rbs(0,1)    =  0.000000000000000E+00 Zbc(0,1)    =  0.000000000000000E+00
Rbc(0,2)    =  2.153736493678703E-02 Zbs(0,2)    =  7.501662133676816E-02 Rbs(0,2)    =  0.000000000000000E+00 Zbc(0,2)    =  0.000000000000000E+00
Rbc(0,3)    = -1.562138963974320E-02 Zbs(0,3)    =  1.387120513836579E-02 Rbs(0,3)    =  0.000000000000000E+00 Zbc(0,3)    =  0.000000000000000E+00
Rbc(0,4)    =  2.560754710143390E-03 Zbs(0,4)    = -5.229796849515394E-03 Rbs(0,4)    =  0.000000000000000E+00 Zbc(0,4)    =  0.000000000000000E+00
Rbc(0,5)    =  3.012855418467466E-04 Zbs(0,5)    = -2.028576061652557E-04 Rbs(0,5)    =  0.000000000000000E+00 Zbc(0,5)    =  0.000000000000000E+00
Rbc(0,6)    = -6.661682495348610E-04 Zbs(0,6)    =  4.167278208885287E-04 Rbs(0,6)    =  0.000000000000000E+00 Zbc(0,6)    =  0.000000000000000E+00
Rbc(0,7)    =  2.411299585436158E-04 Zbs(0,7)    =  1.053846250227302E-04 Rbs(0,7)    =  0.000000000000000E+00 Zbc(0,7)    =  0.000000000000000E+00
Rbc(0,8)    = -4.355168813502804E-05 Zbs(0,8)    = -4.255406393137503E-05 Rbs(0,8)    =  0.000000000000000E+00 Zbc(0,8)    =  0.000000000000000E+00
Rbc(0,9)    =  1.280799530362513E-05 Zbs(0,9)    = -7.001214110899499E-05 Rbs(0,9)    =  0.000000000000000E+00 Zbc(0,9)    =  0.000000000000000E+00
Rbc(0,10)   =  7.054554981072851E-06 Zbs(0,10)   =  3.361702148416995E-05 Rbs(0,10)   =  0.000000000000000E+00 Zbc(0,10)   =  0.000000000000000E+00
Rbc(0,11)   = -1.772234039583128E-05 Zbs(0,11)   =  1.563026464512920E-05 Rbs(0,11)   =  0.000000000000000E+00 Zbc(0,11)   =  0.000000000000000E+00
Rbc(0,12)   =  8.852969544962631E-06 Zbs(0,12)   = -1.533538222669749E-05 Rbs(0,12)   =  0.000000000000000E+00 Zbc(0,12)   =  0.000000000000000E+00
Rbc(0,13)   =  1.838868914837145E-06 Zbs(0,13)   =  5.464039775639624E-08 Rbs(0,13)   =  0.000000000000000E+00 Zbc(0,13)   =  0.000000000000000E+00
Rbc(0,14)   = -3.913528602684541E-06 Zbs(0,14)   =  3.770431023128603E-06 Rbs(0,14)   =  0.000000000000000E+00 Zbc(0,14)   =  0.000000000000000E+00
Rbc(0,15)   =  1.492950910724643E-06 Zbs(0,15)   = -5.597049825772260E-07 Rbs(0,15)   =  0.000000000000000E+00 Zbc(0,15)   =  0.000000000000000E+00
Rbc(0,16)   =  6.159522043139585E-07 Zbs(0,16)   = -3.888037659194830E-07 Rbs(0,16)   =  0.000000000000000E+00 Zbc(0,16)   =  0.000000000000000E+00
Rwc(  0,  0)  =  3.999000000000000E+00 Zws(  0,  0)  =  0.000000000000000E+00 Rws(  0,  0)  =  0.000000000000000E+00 Zwc(  0,  0)  =  0.000000000000000E+00
Rwc(  0,  1)  =  1.800000000000000E+00 Zws(  0,  1)  =  1.800000000000000E+00 Rws(  0,  1)  =  0.000000000000000E+00 Zwc(  0,  1)  =  0.000000000000000E+00
Vns(  0,  0)  =  0.000000000000000E+00 Bns(  0,  0)  =  0.000000000000000E+00 Vnc(  0,  0)  =  1.560452376309866E-15 Bnc(  0,  0)  =  0.000000000000000E+00
Vns(  0,  1)  = -3.794877659599814E-02 Bns(  0,  1)  =  0.000000000000000E+00 Vnc(  0,  1)  = -5.941394574389387E-11 Bnc(  0,  1)  =  0.000000000000000E+00
Vns(  0,  2)  = -3.922424672783899E-02 Bns(  0,  2)  =  0.000000000000000E+00 Vnc(  0,  2)  = -7.258622440325794E-11 Bnc(  0,  2)  =  0.000000000000000E+00
Vns(  0,  3)  =  9.308639008128242E-03 Bns(  0,  3)  =  0.000000000000000E+00 Vnc(  0,  3)  =  2.608500141961564E-10 Bnc(  0,  3)  =  0.000000000000000E+00
Vns(  0,  4)  = -4.028451408453789E-03 Bns(  0,  4)  =  0.000000000000000E+00 Vnc(  0,  4)  =  5.993166721182075E-10 Bnc(  0,  4)  =  0.000000000000000E+00
Vns(  0,  5)  =  3.313589445158859E-03 Bns(  0,  5)  =  0.000000000000000E+00 Vnc(  0,  5)  =  3.576028183065129E-09 Bnc(  0,  5)  =  0.000000000000000E+00
Vns(  0,  6)  = -5.532607439055075E-04 Bns(  0,  6)  =  0.000000000000000E+00 Vnc(  0,  6)  =  3.322463086860999E-09 Bnc(  0,  6)  =  0.000000000000000E+00
Vns(  0,  7)  =  3.428219168629173E-05 Bns(  0,  7)  =  0.000000000000000E+00 Vnc(  0,  7)  =  1.134290954671647E-08 Bnc(  0,  7)  =  0.000000000000000E+00
Vns(  0,  8)  =  8.337673062305028E-05 Bns(  0,  8)  =  0.000000000000000E+00 Vnc(  0,  8)  =  5.287142070379887E-09 Bnc(  0,  8)  =  0.000000000000000E+00
Vns(  0,  9)  = -1.474119492961808E-05 Bns(  0,  9)  =  0.000000000000000E+00 Vnc(  0,  9)  =  9.574223304888734E-09 Bnc(  0,  9)  =  0.000000000000000E+00
Vns(  0, 10)  =  6.218008485843968E-05 Bns(  0, 10)  =  0.000000000000000E+00 Vnc(  0, 10)  =  5.215672560023946E-10 Bnc(  0, 10)  =  0.000000000000000E+00
Vns(  0, 11)  = -5.844833291487409E-05 Bns(  0, 11)  =  0.000000000000000E+00 Vnc(  0, 11)  = -2.256210908750452E-09 Bnc(  0, 11)  =  0.000000000000000E+00
Vns(  0, 12)  =  2.524162293834058E-05 Bns(  0, 12)  =  0.000000000000000E+00 Vnc(  0, 12)  =  8.473406795273612E-10 Bnc(  0, 12)  =  0.000000000000000E+00
Vns(  0, 13)  = -8.506000183907021E-06 Bns(  0, 13)  =  0.000000000000000E+00 Vnc(  0, 13)  = -4.729409425764913E-12 Bnc(  0, 13)  =  0.000000000000000E+00
Vns(  0, 14)  =  2.346231056602112E-06 Bns(  0, 14)  =  0.000000000000000E+00 Vnc(  0, 14)  =  8.465194941860568E-11 Bnc(  0, 14)  =  0.000000000000000E+00
Vns(  0, 15)  =  7.951699281357777E-07 Bns(  0, 15)  =  0.000000000000000E+00 Vnc(  0, 15)  = -7.998112787775051E-11 Bnc(  0, 15)  =  0.000000000000000E+00
Vns(  0, 16)  = -1.745849512900860E-06 Bns(  0, 16)  =  0.000000000000000E+00 Vnc(  0, 16)  = -7.225481421862335E-11 Bnc(  0, 16)  =  0.000000000000000E+00
Vns(  0, 17)  =  1.078840006056463E-06 Bns(  0, 17)  =  0.000000000000000E+00 Vnc(  0, 17)  =  5.592992500901814E-11 Bnc(  0, 17)  =  0.000000000000000E+00
Vns(  0, 18)  = -3.939907951322190E-07 Bns(  0, 18)  =  0.000000000000000E+00 Vnc(  0, 18)  = -5.275497004811282E-12 Bnc(  0, 18)  =  0.000000000000000E+00
Vns(  0, 19)  =  1.217192293403298E-07 Bns(  0, 19)  =  0.000000000000000E+00 Vnc(  0, 19)  = -1.513490248147115E-12 Bnc(  0, 19)  =  0.000000000000000E+00
Vns(  0, 20)  = -1.901651113685301E-08 Bns(  0, 20)  =  0.000000000000000E+00 Vnc(  0, 20)  = -2.880758530145514E-12 Bnc(  0, 20)  =  0.000000000000000E+00
Vns(  0, 21)  = -2.647510301702990E-08 Bns(  0, 21)  =  0.000000000000000E+00 Vnc(  0, 21)  =  5.242906091641308E-14 Bnc(  0, 21)  =  0.000000000000000E+00
Vns(  0, 22)  =  2.780434794984005E-08 Bns(  0, 22)  =  0.000000000000000E+00 Vnc(  0, 22)  =  1.797264411135344E-12 Bnc(  0, 22)  =  0.000000000000000E+00
Vns(  0, 23)  = -1.375025010347072E-08 Bns(  0, 23)  =  0.000000000000000E+00 Vnc(  0, 23)  = -6.732294165505068E-13 Bnc(  0, 23)  =  0.000000000000000E+00
Vns(  0, 24)  =  5.176564574877741E-09 Bns(  0, 24)  =  0.000000000000000E+00 Vnc(  0, 24)  = -1.481945974626548E-14 Bnc(  0, 24)  =  0.000000000000000E+00
Vns(  0, 25)  = -1.825885415747182E-09 Bns(  0, 25)  =  0.000000000000000E+00 Vnc(  0, 25)  = -6.058651251186634E-14 Bnc(  0, 25)  =  0.000000000000000E+00
Vns(  0, 26)  =  1.633239058194226E-10 Bns(  0, 26)  =  0.000000000000000E+00 Vnc(  0, 26)  =  6.357974737647822E-14 Bnc(  0, 26)  =  0.000000000000000E+00
Vns(  0, 27)  =  4.815511727857511E-10 Bns(  0, 27)  =  0.000000000000000E+00 Vnc(  0, 27)  =  2.834478952143090E-14 Bnc(  0, 27)  =  0.000000000000000E+00
Vns(  0, 28)  = -3.830189037666284E-10 Bns(  0, 28)  =  0.000000000000000E+00 Vnc(  0, 28)  = -3.145766694280950E-14 Bnc(  0, 28)  =  0.000000000000000E+00
Vns(  0, 29)  =  1.948323154013029E-10 Bns(  0, 29)  =  0.000000000000000E+00 Vnc(  0, 29)  =  7.877980697803232E-15 Bnc(  0, 29)  =  0.000000000000000E+00
Vns(  0, 30)  = -7.283564805457190E-11 Bns(  0, 30)  =  0.000000000000000E+00 Vnc(  0, 30)  =  4.770337093005714E-16 Bnc(  0, 30)  =  0.000000000000000E+00
Vns(  0, 31)  =  3.566700094412163E-11 Bns(  0, 31)  =  0.000000000000000E+00 Vnc(  0, 31)  =  3.826152854833134E-15 Bnc(  0, 31)  =  0.000000000000000E+00
Vns(  0, 32)  =  9.212623947721024E-12 Bns(  0, 32)  =  0.000000000000000E+00 Vnc(  0, 32)  =  9.446187660745210E-16 Bnc(  0, 32)  =  0.000000000000000E+00
/
&numericlist
 Linitialize =         1
 Lzerovac    =         0
 Ndiscrete   =         2
 Nquad       =        -1
 iMpol       =        -4
 iNtor       =        -4
 Lsparse     =         0
 Lsvdiota    =         1
 imethod     =         3
 iorder      =         2
 iprecon     =         1
 iotatol     =  -1.000000000000000E+00
 Lextrap     =         0
 Mregular    =        -1
/
&locallist
 LBeltrami   =         4
 Linitgues   =         1
/
&globallist
 Lfindzero   =         0
 escale      =   0.000000000000000E+00
 opsilon     =   1.000000000000000E+00
 pcondense   =   4.000000000000000E+00
 epsilon     =   1.000000000000000E-01
 wpoloidal   =   1.000000000000000E+00
 upsilon     =   1.000000000000000E-01
 forcetol    =   1.000000000000000E-12
 c05xmax     =   1.000000000000000E-06
 c05xtol     =   1.000000000000000E-12
 c05factor   =   1.000000000000000E-04
 LreadGF     =         F
 mfreeits    =         1
 gBntol      =   1.000000000000000E-10
 gBnbld      =   0.000000000000000E+00
 vcasingeps  =   1.000000000000000E-12
 vcasingtol  =   1.000000000000000E-06
 vcasingits  =         8
 vcasingper  =         1
/
&diagnosticslist
 odetol      =   1.000000000000000E-07
 nPpts       =       500
 nPtrj       =     -1    -1    -1    -1    -1    -1    -1    -1
 LHevalues   =         F
 LHevectors  =         F
 LHmatrix    =         F
 Lperturbed  =         0
 dpp         =        -1
 dqq         =        -1
 Lcheck      =         1
 Ltiming     =         F
/
&screenlist
 Wpp00aa = T
/
     0     0  4.092972744746912E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  4.091750017338736E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  4.089173197067265E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  4.084479133252834E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  4.076407158020524E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  4.062942717163428E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  4.040456517412868E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
     1     0  1.559901017866803E-01 -2.172697342766468E-01  0.000000000000000E+00  0.000000000000000E+00  3.116675778981305E-01 -4.344431172957540E-01  0.000000000000000E+00  0.000000000000000E+00  4.658600394079777E-01 -6.525791355199720E-01  0.000000000000000E+00  0.000000000000000E+00  6.174197783759314E-01 -8.726049173216391E-01  0.000000000000000E+00  0.000000000000000E+00  7.645993017603614E-01 -1.096198337222746E+00  0.000000000000000E+00  0.000000000000000E+00  9.047516291969493E-01 -1.326323994268314E+00  0.000000000000000E+00  0.000000000000000E+00  1.037511255072165E+00 -1.562461763819107E+00  0.000000000000000E+00  0.000000000000000E+00
     2     0 -1.966349502310396E-04  1.558877324807519E-03  0.000000000000000E+00  0.000000000000000E+00 -7.901487051021386E-04  6.314841387540947E-03  0.000000000000000E+00  0.000000000000000E+00 -1.298365316406752E-03  1.393584296717076E-02  0.000000000000000E+00  0.000000000000000E+00 -1.048013264742250E-03  2.417861889650741E-02  0.000000000000000E+00  0.000000000000000E+00  1.103931328720612E-03  3.688899690205546E-02  0.000000000000000E+00  0.000000000000000E+00  7.145548495686693E-03  5.273845312326476E-02  0.000000000000000E+00  0.000000000000000E+00  2.153736493678703E-02  7.501662133676816E-02  0.000000000000000E+00  0.000000000000000E+00
     3     0 -1.823760611912859E-04  1.863060840595704E-04  0.000000000000000E+00  0.000000000000000E+00 -8.817357460838783E-04  8.873885445533384E-04  0.000000000000000E+00  0.000000000000000E+00 -1.700431385676915E-03  1.411022360899362E-03  0.000000000000000E+00  0.000000000000000E+00 -2.696271510782542E-03  1.664008978374496E-03  0.000000000000000E+00  0.000000000000000E+00 -3.951359567915970E-03  1.629901174884974E-03  0.000000000000000E+00  0.000000000000000E+00 -6.024932354021946E-03  2.164488026646482E-03  0.000000000000000E+00  0.000000000000000E+00 -1.562138963974320E-02  1.387120513836579E-02  0.000000000000000E+00  0.000000000000000E+00
     4     0  4.027030384430831E-05 -5.054864372515776E-05  0.000000000000000E+00  0.000000000000000E+00  1.852882331377017E-04 -2.369570725032600E-04  0.000000000000000E+00  0.000000000000000E+00  4.460716575951569E-04 -5.494730974820281E-04  0.000000000000000E+00  0.000000000000000E+00  8.500518855116362E-04 -9.787907642085662E-04  0.000000000000000E+00  0.000000000000000E+00  1.434036408825270E-03 -1.531782380510813E-03  0.000000000000000E+00  0.000000000000000E+00  2.175629682485349E-03 -2.321729494031129E-03  0.000000000000000E+00  0.000000000000000E+00  2.560754710143390E-03 -5.229796849515394E-03  0.000000000000000E+00  0.000000000000000E+00
     5     0  8.894097102513667E-06 -1.070551725610721E-05  0.000000000000000E+00  0.000000000000000E+00  2.469728445274744E-05 -2.354095273574326E-05  0.000000000000000E+00  0.000000000000000E+00  2.133892924006795E-05 -9.576205671253465E-06  0.000000000000000E+00  0.000000000000000E+00 -1.054857847330186E-05  4.172931622961929E-05  0.000000000000000E+00  0.000000000000000E+00 -8.049062078566358E-05  1.305950669528863E-04  0.000000000000000E+00  0.000000000000000E+00 -1.822805040134370E-04  2.389213096691845E-04  0.000000000000000E+00  0.000000000000000E+00  3.012855418467466E-04 -2.028576061652557E-04  0.000000000000000E+00  0.000000000000000E+00
     6     0 -2.708309471813236E-06  3.354120488618639E-06  0.000000000000000E+00  0.000000000000000E+00 -1.286683858502332E-05  1.503357136671379E-05  0.000000000000000E+00  0.000000000000000E+00 -2.969475193312989E-05  3.177040923844711E-05  0.000000000000000E+00  0.000000000000000E+00 -5.437151618809823E-05  4.730198067328501E-05  0.000000000000000E+00  0.000000000000000E+00 -9.768434862901507E-05  5.391072653865710E-05  0.000000000000000E+00  0.000000000000000E+00 -2.088668089370241E-04  6.002908760408773E-05  0.000000000000000E+00  0.000000000000000E+00 -6.661682495348610E-04  4.167278208885287E-04  0.000000000000000E+00  0.000000000000000E+00
     7     0 -3.894980402114511E-07  4.184507840393577E-07  0.000000000000000E+00  0.000000000000000E+00 -2.425403756881979E-07  1.036686660410108E-08  0.000000000000000E+00  0.000000000000000E+00  2.590265900562383E-06 -3.407610137775769E-06  0.000000000000000E+00  0.000000000000000E+00  1.120601207707956E-05 -9.549221792806400E-06  0.000000000000000E+00  0.000000000000000E+00  3.279687636685329E-05 -1.326749101661353E-05  0.000000000000000E+00  0.000000000000000E+00  9.142507180289504E-05  3.436729511258101E-06  0.000000000000000E+00  0.000000000000000E+00  2.411299585436158E-04  1.053846250227302E-04  0.000000000000000E+00  0.000000000000000E+00
     8     0  2.269607812991237E-07 -2.634091946702981E-07  0.000000000000000E+00  0.000000000000000E+00  1.001633049258848E-06 -1.141544225268696E-06  0.000000000000000E+00  0.000000000000000E+00  1.840608140179453E-06 -1.872922651171651E-06  0.000000000000000E+00  0.000000000000000E+00  1.033749258092494E-06 -1.081875542108147E-06  0.000000000000000E+00  0.000000000000000E+00 -5.719504843876179E-06  3.238731293477604E-06  0.000000000000000E+00  0.000000000000000E+00 -2.832639303177251E-05  1.066489613191852E-05  0.000000000000000E+00  0.000000000000000E+00 -4.355168813502804E-05 -4.255406393137503E-05  0.000000000000000E+00  0.000000000000000E+00
     9     0  3.172258636566504E-08 -3.700113881980589E-08  0.000000000000000E+00  0.000000000000000E+00 -6.635024655075076E-08  7.225969503633539E-08  0.000000000000000E+00  0.000000000000000E+00 -4.340654121757126E-07  4.955014993730635E-07  0.000000000000000E+00  0.000000000000000E+00 -6.717536889691228E-07  5.569218204709642E-07  0.000000000000000E+00  0.000000000000000E+00  9.350296385017240E-07 -2.311929449449482E-06  0.000000000000000E+00  0.000000000000000E+00  8.664153568631331E-06 -1.631986113292113E-05  0.000000000000000E+00  0.000000000000000E+00  1.280799530362513E-05 -7.001214110899499E-05  0.000000000000000E+00  0.000000000000000E+00
    10     0 -2.478458256353884E-08  2.895250605478183E-08  0.000000000000000E+00  0.000000000000000E+00 -9.785177528211289E-08  1.138819194547240E-07  0.000000000000000E+00  0.000000000000000E+00 -1.129931898686006E-07  1.225241552689398E-07  0.000000000000000E+00  0.000000000000000E+00  1.458611085732330E-08  1.893935780114558E-08  0.000000000000000E+00  0.000000000000000E+00  3.533966353413966E-07  3.012151757904196E-07  0.000000000000000E+00  0.000000000000000E+00  1.141101844460620E-06  4.020002300477298E-06  0.000000000000000E+00  0.000000000000000E+00  7.054554981072851E-06  3.361702148416995E-05  0.000000000000000E+00  0.000000000000000E+00
    11     0 -2.206578764852731E-09  2.372789432635515E-09  0.000000000000000E+00  0.000000000000000E+00  1.148927146003187E-08 -1.617091072354190E-08  0.000000000000000E+00  0.000000000000000E+00  6.314528371368071E-08 -6.898027645481700E-08  0.000000000000000E+00  0.000000000000000E+00  3.427262397565998E-08 -3.454860263014715E-08  0.000000000000000E+00  0.000000000000000E+00 -5.137561440883594E-07  4.726316028477491E-07  0.000000000000000E+00  0.000000000000000E+00 -3.494487586342819E-06  2.728362789607470E-06  0.000000000000000E+00  0.000000000000000E+00 -1.772234039583128E-05  1.563026464512920E-05  0.000000000000000E+00  0.000000000000000E+00
    12     0  2.925004055483339E-09 -3.179478303338686E-09  0.000000000000000E+00  0.000000000000000E+00  1.158424502624240E-08 -1.132194311425724E-08  0.000000000000000E+00  0.000000000000000E+00  5.439453899320180E-09 -6.147141900925038E-09  0.000000000000000E+00  0.000000000000000E+00 -2.197042903652287E-09  4.587203187967611E-10  0.000000000000000E+00  0.000000000000000E+00  1.665112859119699E-07 -2.091802726536391E-07  0.000000000000000E+00  0.000000000000000E+00  1.484042059542490E-06 -1.874616686673984E-06  0.000000000000000E+00  0.000000000000000E+00  8.852969544962631E-06 -1.533538222669749E-05  0.000000000000000E+00  0.000000000000000E+00
    13     0  2.040100498504182E-10 -2.170466168827465E-10  0.000000000000000E+00  0.000000000000000E+00 -2.656222300614025E-09  3.119423221928820E-09  0.000000000000000E+00  0.000000000000000E+00 -6.700737441298884E-09  8.573084120732045E-09  0.000000000000000E+00  0.000000000000000E+00  9.920279565815704E-10  3.370494424614381E-09  0.000000000000000E+00  0.000000000000000E+00  5.221166685317585E-08 -1.683329925215095E-08  0.000000000000000E+00  0.000000000000000E+00  1.893520038299817E-07  1.104648337950655E-07  0.000000000000000E+00  0.000000000000000E+00  1.838868914837145E-06  5.464039775639624E-08  0.000000000000000E+00  0.000000000000000E+00
    14     0 -2.385541409232193E-10  4.625036844789103E-10  0.000000000000000E+00  0.000000000000000E+00 -1.151363734492443E-09  7.129946315070713E-10  0.000000000000000E+00  0.000000000000000E+00 -9.793546305524052E-10 -5.461948926040112E-10  0.000000000000000E+00  0.000000000000000E+00 -9.684753128665664E-10 -1.420906985746295E-10  0.000000000000000E+00  0.000000000000000E+00 -3.962625265573144E-08  4.052000535045706E-08  0.000000000000000E+00  0.000000000000000E+00 -3.608274707567829E-07  3.374029901516459E-07  0.000000000000000E+00  0.000000000000000E+00 -3.913528602684541E-06  3.770431023128603E-06  0.000000000000000E+00  0.000000000000000E+00
    15     0 -6.962788879894991E-12  8.395337494774792E-12  0.000000000000000E+00  0.000000000000000E+00  4.879776786031075E-10 -3.373125972511552E-10  0.000000000000000E+00  0.000000000000000E+00  1.210119118000342E-09 -9.013925188997354E-10  0.000000000000000E+00  0.000000000000000E+00  9.603665330942824E-10 -7.354364779554137E-10  0.000000000000000E+00  0.000000000000000E+00  1.507696627547278E-08 -8.165953271185079E-09  0.000000000000000E+00  0.000000000000000E+00  1.842850463489207E-07 -9.989885347162499E-08  0.000000000000000E+00  0.000000000000000E+00  1.492950910724643E-06 -5.597049825772260E-07  0.000000000000000E+00  0.000000000000000E+00
    16     0  5.526160877836174E-10  2.272172848098073E-10  0.000000000000000E+00  0.000000000000000E+00  7.975136012730277E-09  5.838846842019386E-09  0.000000000000000E+00  0.000000000000000E+00  1.525356241529229E-08  1.182470834184808E-08  0.000000000000000E+00  0.000000000000000E+00  2.332867712575383E-08  1.717830380177718E-08  0.000000000000000E+00  0.000000000000000E+00  3.702288577003096E-08  1.854567498303200E-08  0.000000000000000E+00  0.000000000000000E+00  7.914337341705797E-08  1.082780067711484E-08  0.000000000000000E+00  0.000000000000000E+00  6.159522043139585E-07 -3.888037659194830E-07  0.000000000000000E+00  0.000000000000000E+00
zhisong commented 5 years ago

OK. The problem is again coming from the annoying angle definition:

Your input file was made based on the template 'solovev_fb_vmecvol7_final.sp'. This template file was modified from some .end file. In the .end file, the angle convention is defined in SPEC way. But the new Rwc and Zws you put in the input file is the VMEC way.

Changing the sign of Zws will solve the problem instantly.

zhucaoxiang commented 5 years ago

@zhisong You are right. It is working as long as I changed the sign of Zws. Speechless to the angle issue...We need a document to summarize all the angles and conventions used in SPEC.

image

zhisong commented 5 years ago

@zhucaoxiang What did you do to produce the figure in your last post? Was that a converged solution with Lfindzero=2 and free-boundary iteration?

zhucaoxiang commented 5 years ago

@jloizu Yes, that is the figure after following the steps (with flipping the sign of Zws):

  1. run with Lfindzero=0, Lvirtualcasing=inside, gBnbld=0, mfreeits=1 (using Bns=0 as guess)
  2. run with Lfindzero=2, Lvirtualcasing=outside, gBnbld=0, mfreeits=20 (using Bns from previous run)

Please note in the second step, I accidentally set gBnbld=0 (I forgot to change). I will try to put this into an internal loop and revise gBnbld.

jloizu commented 5 years ago

Looks great! I think now we should try a "resolution scan", in order to get something like

|R_axis_spec - R_axis_vmec| versus Nvol

and maybe a similar metric to quantify the "error" in the plasma boundary geometry

zhisong commented 5 years ago

@zhucaoxiang what was the finial bnserr showing in the screen output? I would expect it to be on the order of 1e-9 or 1e-10.

zhucaoxiang commented 5 years ago

@zhisong Yes, at the order of 1E-10.

xspech :    2380.24 : #freeits= 19 ; |f|= 3.81387E-17 ; time=      1.32s ; log|BB|e=-17.44-17.42-17.42-17.50-17.12-17.37-16.84
xspech :            :              ;                                     ; log|II|o=-17.52-17.03-16.62-16.61-16.35-16.25-16.13
xspech :            : 
xspech :    2381.91 : nfreeboundaryiterations =     19 /  00020 ; gBntol = 1.0E-10 ; bnserr = 2.03179E-10 ; bnorml time =      1.67s ;
zhucaoxiang commented 5 years ago

@jloizu @zhisong

Stuart @SRHudson made a convergence study over Nv for benchmarking free-boundary SPEC and VMEC. It is closed to be finished.

In the branch free_bound, I have created a directory at /InputFiles/Verification/FreeBoundVMEC. It contains all the files that Stuart is using for benchmarking. If you have time, please pull the files and check if the results are identical on your end.

jloizu commented 5 years ago

Shall we merge the free-boundary branch with the master branch? I think it is time (it is not a good idea to postpone it unnecessarily, since then merging becomes harder). Or do you think we need to make other modifications before that?

zhisong commented 5 years ago

I agree. Let's push it forward. We need to fine tune this branch before merging. It has two problems, but they are easy to fix.

  1. The initial guess of Bns is constructed using field representation for annulus region. It will not work if we have only one SPEC volume - we need to use the field for toroidal region.
  2. The current code stops working for fixed boundary and Lfindzero=2. @zhucaoxiang Can you make your modification only run for free boundary?
zhucaoxiang commented 5 years ago

@zhisong I have pushed a commit to fix the fixed_boundary Lfindzero=2 bug. Please find a time to test it if it still exists.

I have no clue about the first issue. Is it caused by this part? https://github.com/PrincetonUniversity/SPEC/blob/c2dca08ba1f639081b8dfb19a58d1c19115ac7b4/casing.h#L396-L406

zhucaoxiang commented 5 years ago

Another issue that should be discussed is that if we want to keep the flexibility to let the user decide using the fixed boundary results as initial Bns or not. Right now, the fixed boundary Bns will automatically overid user input Bns. @jloizu @zhisong @SRHudson

abaillod commented 5 years ago

I tested the latest version of the free_bound branch (this includes the changes made by @zhucaoxiang ) on a fixed boundary case (G3V02L1Fi.001.sp in the InputFiles/TestCases directory) and it doesn't work. This is due to this part of the code

  ! run fix_boundary for the first free_boundary iteration
  if (Lfreebound.eq.1 .and. nfreeboundaryiterations.eq.0) then  ! first iteration
     first_free_bound = .true.
     !Mvol = Nvol
     gBnbld_old = gBnbld
     gBnbld = zero
     Lfindzero_old = Lfindzero
     Lfindzero = 0 
     mfreeits_old = mfreeits
     mfreeits = 1
     if (myid.eq.0) write(ounit,'("xspech : ",10X," : First iteration of free boundary calculation : update Bns from plasma.")')
  else
     first_free_bound = .false.
     Mvol = Nvol + Lfreebound
     Lfindzero = Lfindzero_old
     gBnbld = gBnbld_old
     mfreeits = mfreeits_old
  endif

In the case where Lfreebound is not equal to 1, Lfindzero is set to the value of Lfindzero_old - which is not initialized (the same problem stands for mfreeits_old). I solved the problem by moving the lines 

Lfindzero_old = Lfindzero
mfreeits_old = mfreeits

before the beginning of the free boundary iterations. I committed these changes with d2b98ef6536d2b38d22162471255297691b9c721. Please double check / test if this works for you!

jloizu commented 5 years ago

I agree. Let's push it forward. We need to fine tune this branch before merging. It has two problems, but they are easy to fix.

1. The initial guess of Bns is constructed using field representation for annulus region. It will not work if we have only one SPEC volume - we need to use the field for toroidal region.

2. The current code stops working for fixed boundary and Lfindzero=2.
   @zhucaoxiang Can you make your modification only run for free boundary?

I think that point 1 is not really a problem, since we never run free-boundary with one volume. In other words, in free-boundary Mvol=Nvol+1>1. Point 2 has been solved by Antoine Baillod (see previous comment). Please test if this works for you.

We should still discuss about user choice for initialization of Bns.

zhucaoxiang commented 5 years ago

https://github.com/PrincetonUniversity/SPEC/commit/3d469ec9d40a8c42342cb236725c31df7fd04d65#commitcomment-35168832

@jloizu I suppose you have tested this part. I implemented it but moved away without testing.

jloizu commented 5 years ago

3d469ec#commitcomment-35168832

@jloizu I suppose you have tested this part. I implemented it but moved away without testing.

I had tried Antoine's modification, which worked. You then changed it further, and it still works, but actually I don't even know how :) since if Lfreebound=0, then this call is made Lfindzero = Lfindzero_old without having assigned specific values to Lfindzero_old...I think we should change this.

zhucaoxiang commented 5 years ago

@jloizu That's why I made the latest change. Lfindzero = Lfindzero_old will only be executed when Lfreebound == 1 and Lfindzero_old = Lfindzero makes sure that Lfindzero_old has the old value.

jloizu commented 5 years ago

Oh I get it now, sorry, I misread some if/else loops! Makes sense to me.

jloizu commented 5 years ago

Another issue that should be discussed is that if we want to keep the flexibility to let the user decide using the fixed boundary results as initial Bns or not. Right now, the fixed boundary Bns will automatically overid user input Bns. @jloizu @zhisong @SRHudson

Yes, in fact I propose that we introduce some input parameter to decide this. At the very least, this is useful for vacuum field calculations, where we know that Bns=0 from the start; running fixed-boundary first to estimate Bns makes the calculation much longer for nothing. If you think this makes sense, I will introduce such input parameter, something like Linitbns=0,1.

zhucaoxiang commented 5 years ago

@jloizu That could work. I would suggest that the default option is to use fixed-boundary calculations. Unless the user is aware of different options and wants to use customized initialization.

jloizu commented 5 years ago

Done! (see last two commits) Please pull the branch and test, in particular you can test the two new cases I added to InputFiles/TestCases, namely, G3V08L1Fr.001.spand G3V01L0Fr.001.sp. Let me know if it works for you.

zhucaoxiang commented 5 years ago

@jloizu Great, I will check it out later. Another suggestion I might have is actually re-name all the examples with explicit names. Just a comment. We can think about it later.