This is an occasional issue with the inverse lookup. It looks to be an issue with the lookup in the region where theta transitions from 2*pi to zero. If you upload a fixed boundary input file it can be debugged.

The GetBcyl routine in LIBSTELL/Sources/Modules/vmec_utils.fhandles two error cases. If the error flag returns -3 then the point is well outside the equilibrium domain. If the flag returns -1 then there was a failure to converge the inverse lookup to sqrt(ftol)=1E-16. If either of the error codes are thrown a value for S and U may be returned but no computation of B is performed. Now GetBcyl has a parameter fmin_acceptable = 1.E-12_dp which is used to pass points which fail the lookup tollerance but may still be considered 'good'. Using @kudav case I get the following:

So clearly points inside the equilbrium are failing the lookup. These points do return a value for S but not a meaningful Br,Bphi, or Bz.

If I set fmin_acceptable = 1.E-6_dp, then the problem seems to go away:

I'm hesitant to make the change in the code official, but I think it's worth noting this issue. I would suggest that @kudav you modify your sources (line 56 of vmec_utils.f) and recompile.

This issue seems to also exist in W7-X cases using VMEC+VC (this is for an EIM configuration, with the nominal tolerance):

Clearly, the lookup fails here just as in the AUG case. I'm re-running the same case with 1e-6 tolerance to see if the fix is the same.

@kudav Let us know the outcome. The issue is related to the inverse lookup so both runs with VMEC+VC and the -plasma flag (so only VMEC) would have this issue. If changing the tolerance helps then we can make a new pull request.

I could not find any wrong grid points with the decreased tolerance, and the fast ion distribution also seems much less affected by artifacts, so I'd say the fix works but perhaps a more in-depth search for the causes is warranted, since it does have a meaningful impact?

(left: original, right: with decreased tolerance, I also changed the FIDASIM distribution grid, hence the increased resolution)

I'm going to run a scan of the parameter to see what makes the most sense. This is a tollerance on the square of the error for the coordinate inversion. So 1E-6 may roughly translate to 1 mm which is a bit too lax for my tolerance. If 1E-8 or 1E-10 works as well then I'll set the parameter to that and push a new branch.

OK here is the comparrision it seems the magic number is greater than 1.0E-10.

So there are minor but noticeable differences in the fast ion distribution between the tolerance at 1e-6, 1e-8 and 1e-9. No grid points with B_R<2e-5 appear for any of these tolerance, while 30 are found for 1e-12.

There stil exist discontinuities in the derivatives of B_R though (which affect the drifts), and these seem to be relatively independent of the tolerance used:

The discontinuity at the edge is real, at least from the theoretical standpoint. VMEC assumes a surface current at the boundary which guaruntees that it is a flux surface. Essentially this surface current is the representation of the external magnetic field which gives rise to the vacuum component of the field inside the equilbrium. Such an effect is a natural outcome of the VMEC model. One could argue that a free boundary VMEC equilibrium at zero beta which exactly matches a flux surface would have no such feature. The little point in the lower left of the equilbrium is a bit worrying.