Open lazersos opened 2 years ago
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.
To confirm the sharp change in derivatives occurs for s = 0.5/(ns-1)
. So between the half grid and axis.
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.
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.
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.
I recently discovered and diagnosed a bug in the
getBcyl
routine invmec_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 thegetBcyl
routine. I then created a test code which outputs the magnetic field near the magnetic axis region.I created a matlab script to process the file
Here we can see a slight wiggle if you look closely at B. I then made a plot of the derivatives of B.
So clearly there's an issue at the magnetic axis. I suspect this arrises from the following logic (found in
tosuvspace
andflx2cyl
)I think the logic for
jslo
is problematic. The firstIF
statement would prevent theELSE
in the second statement from ever being reached.