Initialization of inner surfaces with Linitialize=1

[zhisong]

For the LHD case I am working on, I found SPEC very difficult to converge (reach force balance) with the Linitialize = 1 option. One problem with this option is that the initial position of the inner surfaces does not change with the input axis position Rac and Zas, in the current version of SPEC.

For example, for two volumes, the inner surface will be initialized only according to the boundary, and has nothing to do with the specified axis.

Detail inspection found the following two lines in preset.h are relevant (line 969 and 970).

!if( im(ii).eq.0 ) then ; psifactor(ii,vvol) = Rscale * tflux(vvol)**half       ! 29 Apr 15;
if( im(ii).eq.0 ) then ; psifactor(ii,vvol) = Rscale * tflux(vvol)**zero       ! 08 Feb 16;

It seems that if we switch the commented lines, the inner surfaces will be initialized according to the specified Rac and Zas. The inner surfaces will be initialized later as

;iRbc(1:mn,vvol) = iRbc(1:mn,0) + ( iRbc(1:mn,lvol) - iRbc(1:mn,0) ) * ( psifactor(1:mn,vvol) / Rscale ) / tflux(lvol)**halfmm(1:mn)
;iZbs(2:mn,vvol) = iZbs(2:mn,0) + ( iZbs(2:mn,lvol) - iZbs(2:mn,0) ) * ( psifactor(2:mn,vvol) / Rscale ) / tflux(lvol)**halfmm(2:mn)

It seems that the second line was introduced on 08 Feb 16. I don't know the reason for these lines.

[jloizu]

@zhisong could you make the change that we discussed on our last zoom meeting? Or I can also do it if you prefer. It would be good to have it committed to GIT so that the master branch properly initializes the surfaces given Rac,Zas guesses.

[jloizu]

Looks like, related to that problem, is the initialization of the coordinates in the innermost volume using Rac information. As it turns out, right now SPEC ignores Rac for constructing the coordinate surfaces in the innermost volume. This does not affect the solution but it can make iota profiles non-sense if the coordinate axis is far from the magnetic axis.

As an example to illustrate this, this the Poincare section obtained for W7-X vacuum (fixed-boundary) calculated with Nvol=1. The magenta dots show the grid points, and the coordinate axis clearly differs from the magnetic axis, even though a good guess for Rac was provided. Also the corresponding iota profile is nonsense (because many surfaces do not go around the coordinate axis).

[jloizu]

Of course, as KAM surfaces move to enforce force-balance, the coordinate axis must be re-evaluated and the initial guess for Rac should become irrelevant. However if Nvol=1or more generally if Lfindzero=0 then no KAM motion occurs and in principle we could keep that guess for the coordinates axis. In other words, IF(Nvol=0 or Lfindzero=0) THEN no need to call rzaxis.

An easy fix is to add an IF statement in packxi.h to avoid such call if Nvol=1 or Lfindzero=0.

The same run now gives this result:

[jloizu]

I am going to create a branch to fix this issue and the one regarding the initialization of the surfaces using Rac explicitly.

[zhisong]

In principle, one will not know the axis precisely until a calculation is made (for example from VMEC). How can one adjust the coordinate axis and avoid the problem in the iota profile when we can't give it a good initial guess? Subroutine rzaxis may not always work fine (the W7X calculation above is a good example).

[jloizu]

Well if the magnetic axis is not known, the best SPEC can do right now is to interpolate inwards from the volume boundary and find a sort of barycenter for the coordinate axis. If, on the other hand, some guess for the axis is known (e.g. from VMEC or from a Biot-Savart calculation in the case of a vacuum), then SPEC should just use this as an initial guess. But we must allow that by, e.g., carrying out the change I proposed.

[zhisong]

What I mean is that, we may come up with a better idea to construct the magnetic axis (hence the iota profile) if there are good flux surfaces. I think it may corresponds to the components of B being zero except for the toroidal component.

Otherwise I agree with you that we should leave the option to prescribe the magnetic axis.

[jloizu]

OK I have created a branch called axisinit in which theses issues are now fixed. I have tested it in W7-X and LHD cases (that were showing these problems) and it works nicely. Also the modifications I made are safe in the sense that they only affect the initialization of the surfaces and innermost coordinates; the pre-conditioning factor remains, so far, the same.

I think we can merge this into master. Please @zhisong try to compile and execute this branch on some test case (e.g. the zero beta multi-volume LHD case you had, but now you can set Linitialize=1 and provide only Rac; or the W7-X test case available on GIT).

[zhisong]

I have tested. The branch is working well on the LHD case I posted above.

[jloizu]

Thanks! I also performed some additional tests on the input test cases on GIT and it all looks good.

I merged the branch into master.

Closing the issue.