Problem with restart in free-boundary

[jloizu]

Usually, when a converged equilibrium is found, one can copy the ".sp.end" output file into the ".sp" input file and set Linitialize=0. Then running the ".sp" input makes SPEC start from a converged equilibrium and thus does not need to iterate on the geometry to reduce the force-imbalance.

I realized that this does not work in free-boundary, unless Lconstraint=-1. By "does not work" I mean that the restarted input has a smaller force imbalance but not zero, and one needs to iterate a couple of times the "copy .end into .sp and run again" procedure in order to (finally) have a ".sp" file that is at equilibrium from the start.

This may be due to how the equilibrium poloidal and toroidal fluxes in the vacuum region are written in the ".end " file. Perhaps not written precisely enough?

[jloizu]

Here is an example of input file (one has to remove the .txt extension) that illustrates what I mean.

l2freeboundaryexample_L0.sp.txt

[jloizu]

As identified by @SRHudson , the problem is that the angle constraint coefficient sweight , which is defined in preset.h, was proportional to tflux in each volume. However, in free-boundary and e.g. Lconstraint=0, the toroidal (and poloidal) flux in the outermost volume, namely the vacuum volume, is modified in order to ensure the constraint of the linking current(s). Thus, when force-balance is obtained, the restart .end file contains a converged geometry but the value of tflux (which is the one initially in the input) has not changed, so that when the code is run again with this new file as input, the part of the force-balance that corresponds to the angle constraints is not converged. For example, one gets initially:

newton : 6.03 : 0 0 ; |f|= 5.43747E-03 ; time= 5.87s ; log|BB|e=-16.44
newton : : ; ; log|II|o= -3.10

[jloizu]

One way to fix this, as proposed by @SRHudson, is to remove the dependency of sweight on tflux,

SALLOCATE( sweight, (1:Mvol), zero )
!sweight(1:Mvol) = upsilon * tflux(1:Mvol)**wpoloidal
do vvol = 1, Mvol ; sweight(vvol) = upsilon * (vvol*one/Nvol)**wpoloidal
enddo

[jloizu]

I created a branch called fixrestart and applied this change. I checked that the code compiles and runs, and that the restart problem has been indeed fixed.

An important caveat is that the angle constraint (and therefore the angle definition itself) is different than before, even if running the same case. That makes comparison with previous runs very hard: the physical solution should not change, but for example the harmonics Rmn,Zmn,etc. will be different.

[jloizu]

One alternative would be to keep the definition of sweight as before and write in the .end file the new values of tflux that have been calculated during the previous run.

[jloizu]

Problem solved via the branch fixrestart. Closing issue.