GetBcyl in vmec_utils.f90 not working near axis

[lazersos]

I recently discovered and diagnosed a bug in the getBcyl routine in vmec_utils.f90. The issue stemmed from a BEAMS3D axisymmetric ITER run where particles near the magnetic axis appeared to not be conserving energy. After much digging into BEAMS3D looking for a bug, I discovered the issue appeared to be stemming from errors in the gradients of the magnetic field. This was determined to be arising from the getBcyl routine. I then created a test code which outputs the magnetic field near the magnetic axis region.

PROGRAM test
	USE read_wout_mod, ONLY: read_wout_file
	USE vmec_utils

	IMPLICIT NONE

    INTEGER          :: ier, i, j
	DOUBLE PRECISION :: r, phi, z, br, bphi, bz, s, u

	INTEGER, PARAMETER :: nr = 256
	INTEGER, PARAMETER :: nz = 256
	DOUBLE PRECISION, PARAMETER :: rmin = 6.1
	DOUBLE PRECISION, PARAMETER :: zmin = -0.4
	DOUBLE PRECISION, PARAMETER :: dr = 0.8
	DOUBLE PRECISION, PARAMETER :: dz = 0.8


	! Read VMEC
	CALL read_wout_file('ITER_15MABURN',ier)

	! Loop
	DO i=1, nr
		DO j = 1, nz
			r = rmin + dr*REAL(i-1)/REAL(nr-1)
			z = zmin + dz*REAL(j-1)/REAL(nz-1)
			phi = 0; ier = 0
			CALL GetBCyl(r,phi,z,br,bphi,bz,s,u,ier)
			WRITE(327,'(2(2X,I5),8(2x,ES20.10),2X,I5)') i,j,r,phi,z,br,bphi,bz,s,u,ier
			CALL FLUSH(6)
		END DO
	END DO

END PROGRAM

I created a matlab script to process the file

%% Make the gradiet plot
vmec_data=read_vmec('wout_ITER_15MABURN.nc');
theta=deg2rad(0:360); zeta=0;
rvmec=cfunct(theta,zeta,vmec_data.rmnc,vmec_data.xm,vmec_data.xn);
zvmec=sfunct(theta,zeta,vmec_data.zmns,vmec_data.xm,vmec_data.xn);
data=importdata('fort.327');
ii=data(:,1);
jj=data(:,2);
nr = max(ii);
nz = max(jj);
r=reshape(data(:,3),[nr nz]);
z=reshape(data(:,5),[nr nz]);
br=reshape(data(:,6),[nr nz]);
bp=reshape(data(:,7),[nr nz]);
bz=reshape(data(:,8),[nr nz]);
b = sqrt(br.*br+bp.*bp+bz.*bz);
s=reshape(data(:,9),[nr nz]);
u=reshape(data(:,10),[nr nz]);
fig=figure('Position',[1 1 1024 768],'Color','white');
subplot(2,3,1);
pixplot(r',z',br'); xlabel('R [m]'); ylabel('Z [m]'); title('B_R'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(2,3,2);
pixplot(r',z',bp'); xlabel('R [m]'); ylabel('Z [m]'); title('B_\phi'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(2,3,3);
pixplot(r',z',bz'); xlabel('R [m]'); ylabel('Z [m]'); title('B_Z'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(2,3,4);
pixplot(r',z',b'); xlabel('R [m]'); ylabel('Z [m]'); title('|B|'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(2,3,5);
pixplot(r',z',s'); xlabel('R [m]'); ylabel('Z [m]'); title('s=(r/a)^2'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(2,3,6);
pixplot(r',z',u'); xlabel('R [m]'); ylabel('Z [m]'); title('u_{VMEC}'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
saveas(fig,'GetBcyl.fig');
saveas(fig,'GetBcyl.png');

Here we can see a slight wiggle if you look closely at B. I then made a plot of the derivatives of B.

%% Now make gradient plots
hr=r(1,2)-r(1,1);
hz=z(2,1)-z(1,1);
[dBdR,dBdZ]=gradient(b,hr,hz);
fig=figure('Position',[1 1 1024*3 768],'Color','white');
subplot(1,3,1);
pixplot(r',z',b'); xlabel('R [m]'); ylabel('Z [m]'); title('|B|'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(1,3,2);
pixplot(r',z',dBdR'); xlabel('R [m]'); ylabel('Z [m]'); title('dB/dr'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
subplot(1,3,3);
pixplot(r',z',dBdZ'); xlabel('R [m]'); ylabel('Z [m]'); title('dZ/dr'); set(gca,'FontSize',14);
hold on; plot(rvmec(1,1),zvmec(1,1),'+r'); plot(rvmec(2,:),zvmec(2,:),'r');
saveas(fig,'diff_GetBcyl.fig');
saveas(fig,'diff_GetBcyl.png');

So clearly there's an issue at the magnetic axis. I suspect this arrises from the following logic (found in tosuvspace and flx2cyl)

!
!     FIND INTERPOLATED s VALUE AND COMPUTE INTERPOLATION WEIGHTS wlo, whi
!     RECALL THAT THESE QUANTITIES ARE ON THE HALF-RADIAL GRID...
!     s-half(j) = (j-1.5)*hs, for j = 2,...ns
!
      hs1 = one/(ns-1)
      jslo = INT(c1p5 + s_in/hs1)
      jshi = jslo+1
      wlo = (hs1*(jshi-c1p5) - s_in)/hs1
      whi = 1 - wlo
      IF (jslo .eq. ns) THEN
!        USE Xhalf(ns+1) = 2*Xhalf(ns) - Xhalf(ns-1) FOR "GHOST" POINT VALUE 1/2hs OUTSIDE EDGE
!        THEN, X = wlo*Xhalf(ns) + whi*Xhalf(ns+1) == Xhalf(ns) + whi*(Xhalf(ns) - Xhalf(ns-1)) 
         jshi = jslo-1
         wlo = 1+whi; whi = -whi
      ELSE IF (jslo .eq. 1) THEN
         jslo = 2
      END IF

!
!     FOR ODD-m MODES X ~ SQRT(s), SO INTERPOLATE Xmn/SQRT(s)
! 
      whi_odd = whi*SQRT(s_in/(hs1*(jshi-c1p5)))
      IF (jslo .ne. 1) THEN
         wlo_odd = wlo*SQRT(s_in/(hs1*(jslo-c1p5)))
      ELSE
         wlo_odd = 0
         whi_odd = SQRT(s_in/(hs1*(jshi-c1p5)))
      END IF

I think the logic for jslo is problematic. The first IF statement would prevent the ELSE in the second statement from ever being reached.

[lazersos]

I created a git repo with the files to re-create this run:
https://github.com/lazersos/getBcyl_test
it's private but just message me to be added to it.

[lazersos]

To confirm the sharp change in derivatives occurs for s = 0.5/(ns-1). So between the half grid and axis.

[lazersos]

More digging, the issue appears to originate from the calculation of Bphi. In the below pictures I've plotted isocontours with 1 equating to ns=2 surface. Thus it's clear that the derivatives of BPHI are incorrect inside of the half grid location.

Adding some more context, the quantities in GetBCyl are calculated like

      Br   = Ru1*bsupu1 + Rv1*bsupv1
      Bphi = R1 *bsupv1
      Bz   = Zu1*bsupu1 + Zv1*bsupv1

Now in this axisymmetric case Rv1=0 and Zv1=0. So from the plots we can see that Ru1, Zu1 must be calculated correctly since they're done on the full grid. Additionally bsupu1 must be correct as well since Br and Bz look OK. R is input to this routine so it's correct as well. Thus the issue must be in the calculation of bsupv1. Oddly it's calculated the same way at bsupu1. However, I think the behavior of B^v should be different near the axis.

[lazersos]

A quick check shows that tosuvspace in read_wout_mod.f90 is not the issue. In the following plots the red line is the inner most gridpoint and we see no evidence of a step,

however if we difference these plots a step becomes apparent at the half grid point,

So it appears the half-grid interpolation in the inner domain in tosuvspace is the source of the issue.

[lazersos]

Update

The issue arises from the algorithm used to integrate the Jacobian. The Jacobian is stored on the half grid because it contains terms like df/ds. However, unlike quantities like Lambda (also on half grid) the way even and odd modes must be treated is different. However, the Jacobian also contains terms like f which have a different even-odd mode behavior at the magnetic axis. Thus interpolation breaks down here for the Jacobian.

The fix is to interpolate R, U and Lambda using the existing algorithms, which work correctly for these quantities. The poloidal and toroidal derivatives are analytic after radial interpolation (dL/du, dL/dv, dR/du, dZ/dv). The evaluate the radial derivatives (dR/ds, dZ/ds) onto the half grid and interpolate with the correct algorithm. The correct interpolation is a work in progress currently.

Once we have this then we can evaluate the Jacobian and B-field anywhere in the VMEC domain correctly. The hope is that we can extend this method to calculate the current density as well. Thus we will only need R, Z, Lambda, dphi/ds, and dchi/ds from the equilibrium to calculate any quantity.