Error in reading gfile ip

[JamieXia]

In line 1788 of openpopcon.py "Ipint = np.abs(np.trapz(y=jtoravg, x=cross_sec_areas))", calculated Ip from this fails to match inputted Ip from gfiles from tokamaker.

[hansec]

@JamieXia how far off is it? The Ip reported by TokaMaker is computed on its grid, which is different from the gEQDSK grid so not surprising if there is some difference.

[matt-pharr]

For reference, here is how the current calculation from the gfile in OpenPOPCON works:

Plasma equilibrium means no net Force on the plasma, so per the momentum equation,

$\rho\left(\frac{ \partial }{ \partial t } +\mathbf{v}\cdot \nabla\right)\mathbf{v} = \mathbf{J}\times \mathbf{B}-\nabla p = 0$

so

$\mathbf{J}\times \mathbf{B}=\nabla p$

or

$\frac{1}{\mu_{0}}\left( \nabla\times \mathbf{B} \right) \times \mathbf{B}=\nabla p$

Following Freidberg, we define the poloidal flux

$\Psi(R,z) = \int_{0}^R\int_{0}^{2\pi}B_{z}(r,z)r d\phi dr=2\pi \int_{0}^{2\pi}B_{z}(r,z) r dr..$

So taking the derivative,

$\frac{1}{2\pi r}\frac{ \partial \Psi }{ \partial r } = B_{z}$

And then, integrating WRT $r$,

$r\frac{ \partial }{ \partial z } \left( \frac{1}{2\pi r}\frac{ \partial \Psi }{ \partial r } \right) = -\frac{ \partial }{ \partial r } \left( rB_{r} \right)$
$\implies -\frac{ \partial }{ \partial z } \left( \frac{1}{2\pi r} \Psi \right) = B_{r}.$

So the poloidal field is simply

$\mathbf{B}_{P} =\frac{1}{r}\nabla \Psi\times \mathbf{e}_{\phi} =\frac{1}{2\pi}\nabla\Psi\times \nabla \phi.$

We then also use Ampere's law again,

$\mu_{0}\mathbf{J}_{T}= \nabla \times \mathbf{B}_P = \nabla\times \left( \frac{1}{r}\nabla \Psi\times \mathbf{e}_{\phi} \right) $
$=\left( \nabla\cdot \left( \frac{1}{r}\nabla\Psi \right) \right) \mathbf{e}_{\phi} $
$=-\frac{1}{r}\left( r\frac{ \partial }{ \partial r } \left( \frac{1}{r}\frac{ \partial \Psi }{ \partial r } \right) +\frac{ \partial ^2\Psi }{ \partial z^2 } \right) \mathbf{e}_{\phi} $
$=-\frac{1}{r}\Delta^*\Psi \mathbf{e}_{\phi} $
$=-\frac{1}{r}\left(-\mu_0 \frac{dp}{d\Psi}r^2 - F \frac{dF}{d\Psi}\right)\mathbf{e}_{\phi}$

where we have defined the stokes operator $\Delta^*$. You can see this implemented here:

OpenPOPCON/src/lib/openpopcon_util.py

Lines 288 to 306 in 6cbe5a6

	zs = np.linspace(gEQDSK['zmid']-0.5*gEQDSK['zdim'], gEQDSK['zmid']+0.5*gEQDSK['zdim'], gEQDSK['psirz'].shape[1])
	rs = np.linspace(gEQDSK['rleft'], gEQDSK['rleft']+gEQDSK['rdim'], gEQDSK['psirz'].shape[0])
	psirz = gEQDSK['psirz']
	psinrz = (psirz-gEQDSK['psimag'])/(gEQDSK['psibry']-gEQDSK['psimag'])
	psin, volgrid, agrid, fs = get_fluxvolumes(gEQDSK, Npsi)
	sqrtpsin = np.sqrt(psin)
	volgrid = np.interp(sqrtpsin,np.sqrt(psin),volgrid)
	ffp = np.asarray(gEQDSK['ffprim'])
	psi = np.linspace(gEQDSK['psimag'],gEQDSK['psibry'],np.shape(ffp)[0])
	fpolrz = np.interp(psirz, psi, gEQDSK['fpol'])
	fprimrz = np.interp(psirz, psi, np.gradient(gEQDSK['fpol'], psi))
	ffprz = np.interp(psirz, psi, gEQDSK['ffprim'])
	pprz = np.interp(psirz, psi, gEQDSK['pprime'])
	jpolrz = 1/rs*fprimrz*np.gradient(psirz, rs,zs)
	gradstarpsi = -(4e-7*np.pi)*rs**2*pprz - 0.5*ffprz
	jtorrz = 1/rs * gradstarpsi

There is one difference that does completely change the current that OpenPOPCON calculates by integration. For your equilibrium @JamieXia (FDC24_TokaMaker_2024-12-05_v16_FS_truncate_true.geqdsk), the flux in the plasma is positive, versus it being negative for gMANTA.

This doesn't make a physics difference as far as I am aware, besides the direction of the field lines, but it does cause the routine issues. Here are the calculated currents each way:

After some poking around, it seems that np.interp is garbage when psi is not strictly increasing. Therefore, for cases with greater flux on the axis than at the boundary, we can just reverse the directions of all interpolated quantities, e.g.

reverse = gEQDSK['psimag'] > gEQDSK['psibry']
ffprz = np.interp(psirz, psi[::1-2*reverse], gEQDSK['ffprim'][::1-2*reverse])

I will open a pull request with this change. The new calculated value of $I_P$ from the supplied gfile is 9.1 MA, which is ~5% off from the true value of 9.6 MA. That is compared to the MANTA equilibrium, where OpenPOPCON calculates 10.002 MA vs a true value of 9.999 MA, only 0.02% off. However, gMANTA has a grid of 129x129 vs the supplied gfile's 65x65, 4x the resolution. @JamieXia could you supply a new version of FDC24_TokaMaker_2024-12-05_v16_FS_truncate_true.geqdsk but with maybe 200x200 resolution? Then this can be tested.