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.

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

Lconstraint=-2, and pyspec higher derivatives of magnetic field #153

Open
wants to merge 29 commits into
base: master
Choose a base branch
from

Conversation

richardnies
Copy link
Collaborator

@richardnies richardnies commented Sep 9, 2021

This branch contains three changes:

  1. the implementation by Elizabeth Paul of Lconstraint=-2 (toroidal flux is adjusted until specified linking current is achieved), extended to work with the new Zernike version of SPEC.

  2. possibility to specify non-zero normal component of the magnetic field in single-volume fixed-boundary calculations (through Vns/Vnc), toggled with inputlist flag Lbdybnzero (by default =.true., i.e. Bsups=0).

  3. py_spec feature: derivatives of the magnetic field with respect to s, zeta and theta can now be obtained.

@ejpaul
Copy link
Collaborator

ejpaul commented Sep 12, 2021

Hi Richard,

Thanks for putting together this PR. A few comments:

  • At the moment, Lconstraint = -2 allows one to both a) impose the normal field in a 1 volume calculation and 2) use the toroidal flux to fix the poloidal linking current. Maybe it would make sense to use a separate input variable for each of these options?
  • It looks like some of the documentation needs to be updated:

    SPEC/ma02aa.f90

    Lines 376 to 385 in ad9f86d

    !> <li> The function \f${\bf f}(\boldsymbol{\mu})\f$, which is computed by mp00ac(), is defined by the input parameter \c Lconstraint:
    !> <ul>
    !> <li> If \c Lconstraint = -1, 0, then \f$\boldsymbol{\mu}\f$ is *not* varied and \c Nxdof=0. </li>
    !> <li> If \c Lconstraint = 1, then \f$\boldsymbol{\mu}\f$ is varied to satisfy the transform constraints;
    !> and \c Nxdof=1 in the simple torus and \c Nxdof=2 in the annular regions.
    !> (Note that in the "simple-torus" region, the enclosed poloidal flux \f$\Delta\psi_p\f$ is not well-defined,
    !> and only \f$\mu=\boldsymbol{\mu}_1\f$ is varied in order to satisfy the transform constraint on the "outer" interface of that volume.) </li>
    !> <li> \todo If \c Lconstraint = 2, then \f$\mu=\boldsymbol{\mu}_1\f$ is varied in order to satisfy the helicity constraint,
    !> and \f$\Delta\psi_p=\boldsymbol{\mu}_2\f$ is *not* varied, and \c Nxdof=1.
    !> (under re-construction)

    SPEC/mp00ac.f90

    Lines 65 to 68 in ad9f86d

    !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& -1 \\
    !> (& 0&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 0 \\
    !> (&{{\,\iota\!\!\!}-}(+1)-\texttt{iota (lvol )}&,& 0&)^T, & \textrm{if }\texttt{Lconstraint} &=& 1 \\
    !> (& ?&,& ?&)^T, & \textrm{if }\texttt{Lconstraint} &=& 2
  • There are a few leftover print statements:

    SPEC/ma02aa.f90

    Line 494 in ad9f86d

    print *,"dtflux(lvol): ", dtflux(lvol)

    SPEC/ma02aa.f90

    Line 556 in ad9f86d

    print *,"End of ma02aa.f90- dtflux(lvol): ", dtflux(lvol)

    SPEC/mp00ac.f90

    Lines 226 to 228 in ad9f86d

    print *,"dtf: ", dtf
    print *,"dpf: ", dpf
    print *,"lmu: ", lmu

    SPEC/mp00ac.f90

    Lines 621 to 624 in ad9f86d

    print *,"curtor (computed): ", dItGpdxtp(0,0,lvol)
    print *,"curpol (computed): ", dItGpdxtp(1,0,lvol)
    print *,"curpol (requested): ", curpol
    print *,"Fdof: ", Fdof(1)

    SPEC/mp00ac.f90

    Line 792 in ad9f86d

    print *,"mu(lvol): ", mu(lvol)
  • At the moment, Lconstraint = -2 only applies to the "plasma" volumes and not the "vacuum" volume. It also appears that we have assumed a single-volume calculation. Can you make sure that the relevant input parameters match these assumptions in global.f90:

    SPEC/global.f90

    Line 1078 in ad9f86d

    FATAL( readin, Lconstraint.lt.-2 .or. Lconstraint.gt.3, illegal Lconstraint )

@richardnies
Copy link
Collaborator Author

Hi Elizabeth,

At the moment, Lconstraint = -2 allows one to both a) impose the normal field in a 1 volume calculation and 2) use the toroidal flux to fix the poloidal linking current. Maybe it would make sense to use a separate input variable for each of these options?

I think the way it is implemented, the normal field can be imposed in a fixed-boundary calculation for any Lconstraint: the change is in preset.f90, where the reading Vns/Vnc is now outside of if(Lfreebound). Am I missing something?

@richardnies richardnies requested a review from jloizu September 29, 2021 17:59
@ejpaul
Copy link
Collaborator

ejpaul commented Sep 30, 2021

Hi Elizabeth,

At the moment, Lconstraint = -2 allows one to both a) impose the normal field in a 1 volume calculation and 2) use the toroidal flux to fix the poloidal linking current. Maybe it would make sense to use a separate input variable for each of these options?

I think the way it is implemented, the normal field can be imposed in a fixed-boundary calculation for any Lconstraint: the change is in preset.f90, where the reading Vns/Vnc is now outside of if(Lfreebound). Am I missing something?

Yes, I agree. I guess at the moment Vns/Vnc are read regardless of the input flags but are defaulted to zero. This seems sensible to me. I would just be careful that there are no other places where Bsups is assumed to be zero. For example,

SPEC/src/curent.f90

Lines 124 to 140 in 0030936

if (Lconstraint .eq. -2) then
call build_vector_potential(lvol, innout, ideriv, 0)
do ii = 1, mn ! loop over Fourier harmonics; 17 May 21;
efmn(ii) = -in(ii)*efmn(ii) ! zeta derivative of Ate
ofmn(ii) = in(ii)*ofmn(ii) ! zeta derivative of Ato
cfmn(ii) = -im(ii)*cfmn(ii) ! theta derivative of Aze
sfmn(ii) = im(ii)*sfmn(ii) ! theta derivative of Azo
enddo ! end of do ii;
call invfft( mn, im(1:mn), in(1:mn), ofmn(1:mn), cfmn(1:mn), sfmn(1:mn), efmn(1:mn), &
Nt, Nz, Bsups(1:Ntz,ideriv), Bsups_2(1:Ntz,ideriv))
Bsups(1:Ntz,ideriv) = Bsups(1:Ntz,ideriv) + Bsups_2(1:ntz,ideriv)
else
Bsups = zero
end if

@richardnies
Copy link
Collaborator Author

I would just be careful that there are no other places where Bsups is assumed to be zero.

Good catch, thanks for pointing that out. Of course, I could just remove "if(Lconstraint .eq. -2)" and always compute Bsups, but explicitly setting it to zero seems cleaner for those cases where we want the boundary to be a flux surface.
Would this be motivation enough to introduce another input variable for the case when the normal BC is non-zero?

Do you (or someone else) perhaps also know how I can check the Doxygen documentation for the branch with the changes I added? I'd like to make sure the formatting worked out.

Thanks,
Richard

@jonathanschilling
Copy link
Collaborator

@richardnies

Do you (or someone else) perhaps also know how I can check the Doxygen documentation for the branch with the changes I added? I'd like to make sure the formatting worked out.

You can run doxygen locally on any branch you like. On GitHub, you only get the documentation from the master branch.
Just type doxygen in the main folder (i.e., where the Doxyfile resides) and you should get a docs folder.
In there, there is html(the website part; examine it by opening index.html in your browser) and latex for the PDF documentation. In the latex folder, you have to do a separate make to get the LaTeX version of the documentation to build. The resulting file should be refman.pdf.

Hope this helps :-)

@ejpaul
Copy link
Collaborator

ejpaul commented Oct 4, 2021

Would this be motivation enough to introduce another input variable for the case when the normal BC is non-zero?

I would be ok either way. Does anyone else have an opinion on this?

@jloizu
Copy link
Collaborator

jloizu commented Oct 11, 2021

I am not going to test (yet) this branch on our machines since I see that you are probably going to make some more changes.

Also, it seems to me that one should really split the "imposing the toroidal flux" option with the "imposiing a non-zero normal field on the boundary" option. The former may be called Lconstraint=-2, and the latter may simply be a flag like Lfreeboundary, e.g. Lboundarybnzero=T by default, and set to false if one wants Vns to be read and imposed at the boundary. Then, if your Lconstraint=-2 option needs to use that always, then a simple "if(Lconstraint=-2) then Lboundarybnzero=F" would suffice. What do you think?

@richardnies
Copy link
Collaborator Author

@richardnies

Do you (or someone else) perhaps also know how I can check the Doxygen documentation for the branch with the changes I added? I'd like to make sure the formatting worked out.

You can run doxygen locally on any branch you like. On GitHub, you only get the documentation from the master branch. Just type doxygen in the main folder (i.e., where the Doxyfile resides) and you should get a docs folder. In there, there is html(the website part; examine it by opening index.html in your browser) and latex for the PDF documentation. In the latex folder, you have to do a separate make to get the LaTeX version of the documentation to build. The resulting file should be refman.pdf.

Hope this helps :-)

Thanks Jonathan! That was very helpful indeed.

I am not going to test (yet) this branch on our machines since I see that you are probably going to make some more changes.

Also, it seems to me that one should really split the "imposing the toroidal flux" option with the "imposiing a non-zero normal field on the boundary" option. The former may be called Lconstraint=-2, and the latter may simply be a flag like Lfreeboundary, e.g. Lboundarybnzero=T by default, and set to false if one wants Vns to be read and imposed at the boundary. Then, if your Lconstraint=-2 option needs to use that always, then a simple "if(Lconstraint=-2) then Lboundarybnzero=F" would suffice. What do you think?

I think you're right, the additional Lboundarybnzero flag is the best way to go. I will implement that now and hopefully we can get this merged then.

…o only allowed with Lconstraint=-2 and Lfreebound=0; and Lconstraint=-2 only allowed with Nvol=1).
@richardnies
Copy link
Collaborator Author

@jloizu @ejpaul I think this PR is now ready.
Thanks,
Richard

@ejpaul
Copy link
Collaborator

ejpaul commented Oct 19, 2021

Sorry I did not see this earlier - I somehow did not receive any notification. I am satisfied with this PR, but I think @jloizu should also give his review.

@jloizu
Copy link
Collaborator

jloizu commented Oct 19, 2021

Sorry I did not review this yet, I will test compilation/execution tomorrow on our machines. Could you add (e.g. on the folder with test cases) an input file where these new options are used?

@richardnies
Copy link
Collaborator Author

@jloizu I have now added the test case in InputFiles/TestCases/G3V01Lm2Fi.info and .sp
Thanks for reviewing!

@richardnies
Copy link
Collaborator Author

It looks like the "py_spec build" check now fails, although I only added two new test case files. Can this be caused by changes I made locally?

@jloizu
Copy link
Collaborator

jloizu commented Oct 21, 2021

Compilation works well on our machines.
However, when testing this case, G2V32L1Fi.001.sp, it does not retrieve the converged equilibrium (force balance ~1e-11) as before. Can you check?
Regarding pyspec build, I have no idea, you need to ask @zhucaoxiang or @mbkumar I guess...

@richardnies
Copy link
Collaborator Author

Compilation works well on our machines.
However, when testing this case, G2V32L1Fi.001.sp, it does not retrieve the converged equilibrium (force balance ~1e-11) as before. Can you check?
Regarding pyspec build, I have no idea, you need to ask @zhucaoxiang or @mbkumar I guess...

Hi @jloizu , how does the lack of force balance manifest itself? I ran the case you quoted with the master and Lconstraint=-2 executables and got very similar outputs, see file attached.

diff_log_G2V32L1Fi.001.txt

@jloizu
Copy link
Collaborator

jloizu commented Oct 22, 2021

For example when I run G3V02L1Fi.001.sp, I get this

xspech :            : version =  3.10
       :  compiled  : date    = Thu Oct 21 14:21:54 CEST 2021 ; 
       :            : srcdir  = /home/loizu/Sandboxes/SPEC ; 
       :            : macros  = src/macros ; 
       :            : fc      = mpif90 ; 
       :            : flags   =  -r8 -DIFORT -O2 -ip -no-prec-div -xHost -fPIC -
 DOPENMP -fopenmp ; 
xspech :            : 
xspech :       0.00 : date=2021/10/22 , 13:32:52 ; machine precision= 1.11E-16 ; vsmall= 1.11E-14 ; small= 1.11E-12 ;
xspech :            : 
xspech :       0.00 : parallelism : ncpu=  1 ; nthreads=  8 ;
rdcmdl :            : 
rdcmdl :       0.00 : ext = G3V02L1Fi.001                                                                                       
xspech :            : 
xspech :       0.00 : begin execution ; calling global:readin ;
readin :            : 
readin :       0.01 : Igeometry=  3 ; Istellsym=  1 ; Lreflect=  0 ; Lbdybnzero= T ;
readin :            : Lfreebound=  0 ; phiedge=  2.000000000000000E+00 ; curtor=  1.038123580200000E-09 ; curpol=  0.000000000000000E+00 ;
readin :            : gamma=  0.000000000000000E+00 ;
readin :            : Nfp=  5 ; Nvol=  2 ; Mvol=  2 ; Mpol=  4 ; Ntor=  4 ;
readin :            : pscale=  1.00000E-03 ; Ladiabatic= 0 ; Lconstraint=  1 ; mupf: tol,its=  1.00E-12 , 128 ;
readin :            : Lrad =  8, 4,
readin :            : 
readin :       0.01 : Linitialize=  0 ;LautoinitBn=  1 ; 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= 3 ;
readin :            : 
readin :       0.01 : LBeltrami= 4 ; Linitgues= 1 ; Lmatsolver= 3 ; LGMRESprec= 1 ; NiterGMRES= 200 ; epsGMRES=  1.00000E-14 ; epsILU=  1.00000E-12 ;
readin :            : 
readin :       0.01 : Lfindzero= 2 ;
readin :            : escale=  0.00000E+00 ; opsilon=  1.00000E+00 ; pcondense=  4.000 ; epsilon=  1.00000E+00 ; wpoloidal= 1.0000 ; upsilon=  1.00000E+00 ;
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-06 ; gBnbld=  6.66000E-01 ;
readin :            : vcasingeps=  1.00000E-12 ; vcasingtol=  1.00000E-08 ; vcasingits=     8 ; vcasingper=     1 ;
readin :            : 
readin :       0.01 : odetol=  1.00E-07 ; nPpts=     0 ;
readin :            : LHevalues= F ; LHevectors= F ; LHmatrix= F ; Lperturbed= 0 ; dpp= -1 ; dqq= -1 ; dRZ=  1.00000000E-05 ; Lcheck=  0 ; Ltiming= F ;
readin :            : 
readin :            : myid=  0 ; Rscale= 1.000000000000000E+01 ;
preset :            : myid=  0 ; Mrad=  8 : Lrad=  8,  4,
preset :       0.01 : LBsequad= F , LBnewton= F , LBlinear= T ;
preset :            : 
preset :       0.01 : Nquad=  -1 ; mn=   41 ; NGdof=    81 ; NAdof=   361,   370,
preset :            : 
preset :       0.01 : Nt=    32 ; Nz=    32 ; Ntz=     1024 ;
mp00ac :       0.10 : myid=  0 ; lvol=  1 ; ideriv= 0 ; GMRES ierr=  -3 ;      solver internal break down ;  
mp00ac :       0.15 : myid=  0 ; lvol=  1 ; ideriv= 0 ; GMRES ierr=  -3 ;      solver internal break down ;  
mp00ac :       0.15 : myid=  0 ; lvol=  1 ; ideriv= 1 ; GMRES ierr=  -3 ;      solver internal break down ;  
mp00ac :       0.15 : myid=  0 ; lvol=  1 ; ideriv= 0 ; GMRES ierr=  -3 ;      solver internal break down ;  
tr00ab :       0.15 ; myid=  0 ; lvol=  1 ; innout= 1 ; jderiv= 2 ; idgesvx= 1 ; time=    0.0001 ; singular ;       
tr00ab :      fatal : myid=  0 ; idgesvx.ne.0 ; failed to construct straight-fieldline angle using dgesvx ;

while with an older version I get this

xspech :            : version =  3.01
       :  compiled  : date    = Fri Jun 25 12:03:54 CEST 2021 ; 
       :            : srcdir  = /home/loizu/Sandboxes/SPEC ; 
       :            : macros  = macros ; 
       :            : fc      = mpif90 ; 
       :            : flags   =  -r8 -O2 -ip -no-prec-div -xHost -fPIC -DOPENMP 
 -fopenmp ; 
xspech :            : 
xspech :       0.00 : begin execution ; ncpu=  1 ; omp_threads=  8; calling global:readin ;
readin :            : 
readin :       0.00 : date=2021/10/22 , 13:32:30 ; machine precision= 1.11E-16 ; vsmall= 1.11E-14 ; small= 1.11E-12 ;
readin :            : 
readin :       0.00 : ext = G3V02L1Fi.001                                                                                       
readin :            : 
readin :            : 
readin :       0.00 : Igeometry=  3 ; Istellsym=  1 ; Lreflect=  0 ;
readin :            : Lfreebound=  0 ; phiedge=  2.000000000000000E+00 ; curtor=  1.038123580200000E-09 ; curpol=  0.000000000000000E+00 ;
readin :            : gamma=  0.000000000000000E+00 ;
readin :            : Nfp=  5 ; Nvol=  2 ; Mvol=  2 ; Mpol=  4 ; Ntor=  4 ;
readin :            : pscale=  1.00000E-03 ; Ladiabatic= 0 ; Lconstraint=  1 ; mupf: tol,its=  1.00E-12 , 128 ;
readin :            : Lrad =  8, 4,
readin :            : 
readin :       0.00 : Linitialize=  0 ;LautoinitBn=  1 ; 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= 3 ;
readin :            : 
readin :       0.00 : LBeltrami= 4 ; Linitgues= 1 ; Lmatsolver= 3 ; LGMRESprec= 1 ; NiterGMRES= 200 ; epsGMRES=  1.00000E-14 ; epsILU=  1.00000E-12 ;
readin :            : 
readin :       0.00 : Lfindzero= 2 ;
readin :            : escale=  0.00000E+00 ; opsilon=  1.00000E+00 ; pcondense=  4.000 ; epsilon=  1.00000E+00 ; wpoloidal= 1.0000 ; upsilon=  1.00000E+00 ;
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-06 ; gBnbld=  6.66000E-01 ;
readin :            : vcasingeps=  1.00000E-12 ; vcasingtol=  1.00000E-08 ; vcasingits=     8 ; vcasingper=     1 ;
readin :            : dxdesc=  1.00000E-03 ; ftoldesc=  1.00000E-10 ; maxitdesc= 15000 ; Lwritedesc= 1 ; nwritedesc= 100 ; Manderson= 1 ;
readin :            : 
readin :       0.00 : odetol=  1.00E-07 ; nPpts=     0 ;
readin :            : LHevalues= F ; LHevectors= F ; LHmatrix= F ; Lperturbed= 0 ; dpp= -1 ; dqq= -1 ; dRZ=  1.00000000E-05 ; Lcheck=  0 ; Ltiming= F ;
readin :            : 
readin :            : myid=  0 ; Rscale= 1.000000000000000E+01 ;
preset :            : myid=  0 ; Mrad=  8 : Lrad=  8,  4,
preset :       0.00 : LBsequad= F , LBnewton= F , LBlinear= T ;
preset :            : 
preset :       0.00 : Nquad=  -1 ; mn=   41 ; NGdof=    81 ; NAdof=   361,   370,
preset :            : 
preset :       0.00 : Nt=    32 ; Nz=    32 ; Ntz=     1024 ;
newton :       0.07 :         0  0 ; |f|= 1.48437E-02 ; time=      0.06s ; log|BB|e= -4.66
newton :            :              ;                                     ; log|II|o= -2.42
fcn2   :       0.57 :         1  1 ; |f|= 1.48437E-02 ; time=      0.56s ; log|BB|e= -4.66
fcn2   :            :              ;                                     ; log|II|o= -2.42
fcn2   :       0.61 :         2  1 ; |f|= 1.46506E-02 ; time=      0.04s ; log|BB|e= -4.74
fcn2   :            :              ;                                     ; log|II|o= -2.42
fcn2   :       0.65 :         3  1 ; |f|= 1.42655E-02 ; time=      0.04s ; log|BB|e= -4.96
fcn2   :            :              ;                                     ; log|II|o= -2.44
fcn2   :       0.69 :         4  1 ; |f|= 1.35002E-02 ; time=      0.04s ; log|BB|e= -4.65
fcn2   :            :              ;                                     ; log|II|o= -2.47
fcn2   :       0.73 :         5  1 ; |f|= 1.19909E-02 ; time=      0.04s ; log|BB|e= -4.15
fcn2   :            :              ;                                     ; log|II|o= -2.53
fcn2   :       0.77 :         6  1 ; |f|= 9.08607E-03 ; time=      0.04s ; log|BB|e= -3.81
fcn2   :            :              ;                                     ; log|II|o= -2.64
fcn2   :       0.81 :         7  1 ; |f|= 4.32438E-03 ; time=      0.04s ; log|BB|e= -3.57
fcn2   :            :              ;                                     ; log|II|o= -2.85
fcn2   :       0.85 :         8  1 ; |f|= 1.21257E-03 ; time=      0.04s ; log|BB|e= -4.24
fcn2   :            :              ;                                     ; log|II|o= -3.25
fcn2   :       0.91 :         9  1 ; |f|= 2.55751E-03 ; time=      0.06s ; log|BB|e= -4.65
fcn2   :            :              ;                                     ; log|II|o= -2.91
fcn2   :       0.96 :        10  1 ; |f|= 7.31985E-04 ; time=      0.05s ; log|BB|e= -4.45
fcn2   :            :              ;                                     ; log|II|o= -3.41
fcn2   :       1.02 :        11  1 ; |f|= 1.44780E-03 ; time=      0.06s ; log|BB|e= -4.47
fcn2   :            :              ;                                     ; log|II|o= -3.03
fcn2   :       1.07 :        12  1 ; |f|= 5.42314E-04 ; time=      0.05s ; log|BB|e= -4.66
fcn2   :            :              ;                                     ; log|II|o= -3.48
fcn2   :       1.12 :        13  1 ; |f|= 4.39564E-04 ; time=      0.05s ; log|BB|e= -4.79
fcn2   :            :              ;                                     ; log|II|o= -3.53
fcn2   :       1.27 :        14  1 ; |f|= 8.33093E-05 ; time=      0.16s ; log|BB|e= -5.29
fcn2   :            :              ;                                     ; log|II|o= -4.29
fcn2   :       1.35 :        15  1 ; |f|= 6.10494E-05 ; time=      0.07s ; log|BB|e= -5.52
fcn2   :            :              ;                                     ; log|II|o= -4.39
fcn2   :       1.39 :        16  1 ; |f|= 1.33442E-05 ; time=      0.05s ; log|BB|e= -6.18
fcn2   :            :              ;                                     ; log|II|o= -5.06
fcn2   :       1.44 :        17  1 ; |f|= 7.78188E-06 ; time=      0.04s ; log|BB|e= -6.70
fcn2   :            :              ;                                     ; log|II|o= -5.30
fcn2   :       1.48 :        18  1 ; |f|= 4.16587E-06 ; time=      0.04s ; log|BB|e= -6.63
fcn2   :            :              ;                                     ; log|II|o= -5.69
fcn2   :       1.51 :        19  1 ; |f|= 2.58076E-06 ; time=      0.04s ; log|BB|e= -6.85
fcn2   :            :              ;                                     ; log|II|o= -5.89
fcn2   :       1.55 :        20  1 ; |f|= 8.29873E-08 ; time=      0.04s ; log|BB|e= -8.33
fcn2   :            :              ;                                     ; log|II|o= -7.30
fcn2   :       1.58 :        21  1 ; |f|= 1.95748E-08 ; time=      0.03s ; log|BB|e= -8.95
fcn2   :            :              ;                                     ; log|II|o= -7.94
fcn2   :       1.61 :        22  1 ; |f|= 5.43977E-10 ; time=      0.03s ; log|BB|e=-10.51
fcn2   :            :              ;                                     ; log|II|o= -9.43
fcn2   :       1.64 :        23  1 ; |f|= 7.17071E-11 ; time=      0.03s ; log|BB|e=-11.38
fcn2   :            :              ;                                     ; log|II|o=-10.27
fcn2   :       1.66 :        24  1 ; |f|= 9.22314E-12 ; time=      0.02s ; log|BB|e=-12.28
fcn2   :            :              ;                                     ; log|II|o=-11.09
fcn2   :       1.68 :        25  1 ; |f|= 3.25930E-12 ; time=      0.02s ; log|BB|e=-12.72
fcn2   :            :              ;                                     ; log|II|o=-11.55
fcn2   :       1.71 :        26  1 ; |f|= 1.34467E-13 ; time=      0.02s ; log|BB|e=-14.03
fcn2   :            :              ;                                     ; log|II|o=-12.93
fcn2   :       1.72 :        27  1 ; |f|= 2.85433E-14 ; time=      0.02s ; log|BB|e=-14.80
fcn2   :            :              ;                                     ; log|II|o=-13.64
fcn2   :       1.74 :        28  1 ; |f|= 9.64711E-15 ; time=      0.02s ; log|BB|e=-15.12
fcn2   :            :              ;                                     ; log|II|o=-14.06
fcn2   :       1.75 :        29  1 ; |f|= 2.39385E-15 ; time=      0.01s ; log|BB|e=-15.61
fcn2   :            :              ;                                     ; log|II|o=-14.75
newton :            :
newton :       1.75 : finished ; success        ; ic05p*f= 1 ; its=     29 ,   1 ;
xspech :       1.76 : #freeits=  0 ; |f|= 2.39385E-15 ; time=      0.01s ; log|BB|e=-15.61
xspech :            :              ;                                     ; log|II|o=-14.75
xspech :            :
xspech :       1.79 : myid=  0 : time=    0.03m =   0.00h =  0.00d ;
ending :            : 
ending :       1.79 : myid=  0 ; completion ; time=      1.79s =     0.03m =   0.00h =  0.00d ; date= 2021/10/22 ; time= 13:32:32 ; ext = G3V02L1Fi.001                                               
ending :            : 

@richardnies
Copy link
Collaborator Author

Hi @jloizu ,
are you able to run that same case G3V02L1Fi.001.sp with a new xspec from the master branch? I can't run this case (on master or Lconstraint_-2 branches) due to the error "volume : fatal : myid= 0 ; vflag.eq.0 .and. vvolume(lvol).lt.small ; volume cannot be zero or negative", which I've encountered for other cases before.

@jloizu
Copy link
Collaborator

jloizu commented Oct 22, 2021

Indeed there seems to be something already screwed up in the master branch...this is not good.

I don't even understand how the previous versions passed the CI tests, since I just tried running some of them with the master branch and the code does not find equilibria...We should open an issue on this!

So my guess is that your changes and branch are all fine, but there is some other more general issue...

@jonathanschilling
Copy link
Collaborator

@jloizu @richardnies The problem volume : fatal : myid= 0 ; vflag.eq.0 .and. vvolume(lvol).lt.small ; volume cannot be zero or negative arises on my machine if I build SPEC using the CMake setup. When built using the traditional Makefile, it works.
The problem seems to be that in the current CMake setup the NSCREENLIST macro is not expanded.
Thus in G3V02L1Fi.001.sp reading the screenlist fails since Wpp00aa is not member of the namelist screenlist.

@mbkumar Could you please fix the expansion of NSCREENLIST in the CMake setup?

I added a few checks on the master branch in global.f90:read_inputlists_from_file() to help debug further problems like this.

@jonathanschilling
Copy link
Collaborator

@jloizu

I don't even understand how the previous versions passed the CI tests, since I just tried running some of them with the master branch and the code does not find equilibria...We should open an issue on this!

I think this could be because in the CI on GitLab Actions we still build SPEC using the Makefile setup, but @richardnies probably uses the CMake setup (as preferred these days)?

@mbkumar
Copy link
Collaborator

mbkumar commented Oct 25, 2021 via email

@richardnies
Copy link
Collaborator Author

richardnies commented Oct 26, 2021

@jonathanschilling

@jloizu

I don't even understand how the previous versions passed the CI tests, since I just tried running some of them with the master branch and the code does not find equilibria...We should open an issue on this!

I think this could be because in the CI on GitLab Actions we still build SPEC using the Makefile setup, but @richardnies probably uses the CMake setup (as preferred these days)?

Actually I use the Makefile setup as well. With gfortran, if that helps.

EDIT: I just saw you updated the master branch, the test runs fine now. Thanks!

@mbkumar
Copy link
Collaborator

mbkumar commented Oct 26, 2021 via email

@jonathanschilling
Copy link
Collaborator

@ejpaul @richardnies

EDIT: I just saw you updated the master branch, the test runs fine now. Thanks!

I re-ran the suggested test case just now and they still seem to work. I guess we can merge this then?

@richardnies
Copy link
Collaborator Author

@jonathanschilling Sorry for not getting back to you earlier. One change that happened since the code changes were reviewed is that I implemented FFT's to obtain the B-field in py_spec. Should that be reviewed by someone?

Apart from that, it is ready to merge on my side.

@jonathanschilling
Copy link
Collaborator

@richardnies No worries :-)
I like the implementation of FFTs and I would be happy to test this.
Could you maybe include an example script in the repository, which shows how to use these new routines along with some plots/... of expected results?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

5 participants