Read initial guess with python wrappers

[abaillod]

Hi,

I want to read and modify the initial guess from an input file using the python wrappers. The goal is to implement a new feature in simsopt, where the latest converged equilibrium would be used as initial guess for the next equilibrium calculation. However, it does not seem to work - after running spec.allglobal.read_inputlists_from_file(), the arrays spec.allglobal.allrzrz, spec.allglobal.mmrzrz and spec.allglobal.nnrzrz are empty.

As a minimal working example, consider the input file test.txt (I had to change the file extension to attach it, please change it to .sp when you run it). Running the following python script:

import spec.spec_f90wrapped as spec
from mpi4py import MPI 
import numpy as np

spec.inputlist.initialize_inputs()
spec.allglobal.ext = 'test'
spec.allglobal.read_inputlists_from_file()
spec.allglobal.check_inputs()

print('------------------------------------------------------------')
print(' CHECKING INITIAL GUESS: ')
print('allrzrz = ', spec.allglobal.allrzrz)
print('mmrzrz = ', spec.allglobal.mmrzrz)
print('nnrzrz = ', spec.allglobal.nnrzrz)

print('------------------------------------------------------------')
print(' WRITING INITIAL GUESS: ')
spec.allglobal.allrzrz = np.zeros((4,2,25))

generates the following output and error message:

(spec_env) abaillod@spcpc602:~/development/sandbox/Optimization/test_spec_initial_guess> python script.py 
Invalid line in screenlist:      2     2  1.323348147776832E-03  1.164151596410146E-03  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
readin :            : 
readin : ********** : Igeometry=  3 ; Istellsym=  1 ; Lreflect=  0 ;
readin :            : Lfreebound=  0 ; phiedge=  1.000000000000000E+00 ; curtor=  0.000000000000000E+00 ; curpol=  2.143330000000000E+01 ;
readin :            : gamma=  0.000000000000000E+00 ;
readin :            : Nfp=  5 ; Nvol=  2 ; Mvol=  2 ; Mpol=  2 ; Ntor=  2 ;
readin :            : pscale=  1.84084E-01 ; Ladiabatic= 0 ; Lconstraint=  3 ; mupf: tol,its=  1.00E-12 , 128 ;
readin :            : Lrad =  6, 6,
readin :            : 
readin : ********** : Linitialize=  0 ;LautoinitBn=  0 ; Lzerovac= 0 ; Ndiscrete= 2 ;
readin :            : Nquad=  -1 ; iMpol=  -4 ; iNtor=  -4 ;
readin :            : Lsparse= 0 ; Lsvdiota= 0 ; imethod= 3 ; iorder= 2 ; iprecon= 1 ; iotatol= -1.00000E+00 ;
readin :            : Lextrap= 0 ; Mregular= -1 ; Lrzaxis= 1 ; Ntoraxis= 2 ;
readin :            : 
readin : ********** : LBeltrami= 4 ; Linitgues= 1 ; Lmatsolver= 3 ; LGMRESprec= 1 ; NiterGMRES= 200 ; epsGMRES=  1.00000E-14 ; epsILU=  1.00000E-12 ;
readin :            : 
readin : ********** : Lfindzero= 2 ;
readin :            : escale=  0.00000E+00 ; opsilon=  1.00000E+00 ; pcondense=  4.000 ; epsilon=  1.00000E-01 ; wpoloidal= 1.0000 ; upsilon=  1.00000E-01 ;
readin :            : forcetol=  1.00000E-12 ; c05xmax=  1.00000E-06 ; c05xtol=  1.00000E-12 ; c05factor=  1.00000E-04 ; LreadGF= F ; 
readin :            : mfreeits=   0 ; gBntol=  1.00000E-05 ; gBnbld=  6.00000E-01 ;
readin :            : vcasingeps=  0.00000E+00 ; vcasingtol=  1.00000E-06 ; vcasingits=     0 ; vcasingper=     1 ;
readin :            : 
readin : ********** : odetol=  1.00E-07 ; nPpts=   100 ;
readin :            : LHevalues= F ; LHevectors= F ; LHmatrix= F ; Lperturbed= 0 ; dpp= -1 ; dqq= -1 ; dRZ=  1.00000000E-05 ; Lcheck=  0 ; Ltiming= F ;
readin :            : 
------------------------------------------------------------
 CHECKING INITIAL GUESS: 
allrzrz =  []
mmrzrz =  []
nnrzrz =  []
------------------------------------------------------------
 WRITING INITIAL GUESS: 
Traceback (most recent call last):
  File "/home/abaillod/development/sandbox/Optimization/test_spec_initial_guess/script.py", line 23, in <module>
    spec.allglobal.allrzrz = np.zeros((4,2,25))
  File "/misc/anaconda3/envs/spec_env/lib/python3.10/site-packages/spec/spec_f90wrapped.py", line 5511, in allrzrz
    self.allrzrz[...] = allrzrz
ValueError: could not broadcast input array from shape (4,2,25) into shape (4,2,0)

You see that (i) the arrays which should contain the initial guess harmonics are empty and (ii) trying to modify them / change their size generates an error.

In summary,

• How do I read an initial guess using the python wrappers ? And why what am I doing in the example above does not work ?
• How do I modify the initial guess / change the number of modes (size of the array) ?

I would appreciate any help ! Thank you

[jonathanschilling]

Hi @abaillod ,

the problem is detected in the first line already:

Invalid line in screenlist: 2 2 1.323348147776832E-03 1.164151596410146E-03 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00

SPEC expects every namelist to be present in the input file, even if it is empty. screenlist is missing and thus the input reader tries to interpret the first line of the interface geometry initial guess as screenlist, which fails with above error message. Thus, you have to put the following into the input file between diagnosticslist and the initial guess for the interface geometry:

&screenlist
/

Then it should work... :-)

[abaillod]

Well, I guess I was not paying attention enough... Thanks!

[abaillod]

Hello again,

I re-open this issue because I think there is a type conversion error in the python wrappers when reading the initial guess. Using the same input file as mentionned in my first message, and running the following script:

import spec.spec_f90wrapped as spec
from mpi4py import MPI 
import numpy as np

spec.inputlist.initialize_inputs()
spec.allglobal.ext = 'test'
spec.allglobal.read_inputlists_from_file()
spec.allglobal.check_inputs()

print('------------------------------------------------------------')
print(' CHECKING INITIAL GUESS: ')
print('allrzrz = ', spec.allglobal.allrzrz)
print('mmrzrz = ', spec.allglobal.mmrzrz)
print('nnrzrz = ', spec.allglobal.nnrzrz)

leads to:

... some outputs ...
------------------------------------------------------------
 CHECKING INITIAL GUESS: 
allrzrz =  [[[ 1.02713088e+01  5.15105071e-02  1.10731422e-04  1.19642014e-04
   -9.06594710e-04  4.64232728e-01  1.02925926e-01 -1.41000555e-03
   -1.05129748e-05 -2.00360058e-04  2.60667759e-02  5.36078239e-03
    1.32334815e-03]
  [ 1.00000000e+01  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  1.00000000e+00  2.50000000e-01  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00]]

 [[ 0.00000000e+00  4.82264437e-02  6.74426935e-05 -1.14611422e-04
    1.72699027e-03 -5.45037848e-01  1.15678588e-01 -1.57134537e-03
   -1.16772309e-05  3.34396978e-04 -2.20819071e-02 -4.55664612e-03
    1.16415160e-03]
  [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00 -1.00000000e+00  2.50000000e-01  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00]]

 [[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00]
  [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00]]

 [[ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00]
  [ 0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00  0.00000000e+00  0.00000000e+00  0.00000000e+00
    0.00000000e+00]]]
mmrzrz =  [             0     4294967296     4294967297     4294967297
     8589934594     8589934594 47343424503810            657
              0              0              0 93945306052672
             33]
nnrzrz =  [    4294967296    -8589934590     4294967295     8589934593
             -2     4294967296 47343424503810            657
              0              0              0 93945306052672
             33]

You see that the Rbc, Rbs, Zbc, Zbs have the right harmonics while the values of the mode numbers m and n are wrong. Printing instead print('mmrzrz = ', spec.allglobal.mmrzrz.astype('int16')) gives

mmrzrz =  [     0      0      1      1      2      2      2    657      0      0      0 -10176     33]

which is closer to what we should get (i.e. 0 0 1 1 1 1 1 2 2 2 2 2 ) but still wrong. This makes me think it is probably a type conversion error in the python wrappers. Could anyone help me debug this ?

Thanks !

[zhucaoxiang]

@abaillod Try the branch pywrapper_fix. I have pushed some changes. It is now fixed on my end. But we should be careful that it does not import new bugs.

In [5]: print('mmrzrz = ', spec.allglobal.mmrzrz)
mmrzrz =  [0 0 0 1 1 1 1 1 2 2 2 2 2]

[abaillod]

Thank you @zhucaoxiang, it seems to fix the problem on my side as well.

Could we merge these fixes with the master branch ? I would like to use these in a future simsopt PR.