Helicity must be set in gpec namelist if not using coils

[logan-nc]

The GPEC OMFIT workflow manager previously automated the setting of ip_direction and bt_direction in the coil namelist, but @ysm1887 and @jmhanson found that the user needed to manually set the directions in the gpec namelist.

This is a record of that discussion, so future users/developers can reference the documented issue.

From @jmhanson:

I do have one question about the plasma helicity. 152789 has Ip > 0 and Bt < 0, so I would call the helicity "left-handed". In gpec_control_output_n1.nc, the helicity attribute is +1. I also see that the response mode contours (pages 4-5) go the opposite direction I would expect on the LFS midplane. I think GPEC is getting the helicity from the geqdsk file. Could this information be getting lost or reversed somewhere?

@logan-nc:

GPEC takes the ip_direction and bt_direction as manual inputs in either the COIL_CONTROL namelist or the GPEC_INPUT namelist. OMFIT does indeed attempt to get this information from the gfile as a post-process of the GUI for loading equilibrium, but it only does so through the coil namelist (i.e. the set_coil_equilibrium script). Is it possible you ran this (a) without using the GUI and without calling this script or (b) without using coils but without adding the helicity information to the gpec namelist?

@logan-nc then added this logic to double-check the helicity settings at runtime in the OMFIT run_exes.py script

        # check helicity logic has been properly set
        btd = root['INPUTS'][f'{gexeu}'][f'{gexeu}_INPUT'].get('bt_direction', 'positive')
        ipd = root['INPUTS'][f'{gexeu}'][f'{gexeu}_INPUT'].get('ip_direction', 'positive')
        btip_nml = gexeu
        if root['INPUTS'][gexeu][gexeu + '_INPUT']['coil_flag']:  # coil namelist
            btd = root['INPUTS']['COIL'][t]['COIL_CONTROL'].get('bt_direction', btd)
            ipd = root['INPUTS']['COIL'][t]['COIL_CONTROL'].get('ip_direction', ipd)
            btip_nml = 'COIL'
            if btd != root['INPUTS'][f'{gexeu}'][f'{gexeu}_INPUT'].get('bt_direction', btd) or ipd != root['INPUTS'][f'{gexeu}'][
                f'{gexeu}_INPUT'
            ].get('ip_direction', ipd):
                raise OMFITexception(f"{gexeu} ip & bt directions do not match COIL namelist")
        if root['INPUTS']['EQUIL']['EQUIL_CONTROL']['eq_type'] == 'efit':
            eq = root['INPUTS']['EQUILIBRIUM'][t]
            if (eq['BCENTR'] > 0 and btd != 'positive') or (eq['BCENTR'] < 0 and btd != 'negative'):
                raise OMFITexception(f"{btip_nml} bt_direction does not match g-file")
            if (eq['CURRENT'] > 0 and ipd != 'positive') or (eq['CURRENT'] < 0 and ipd != 'negative'):
                raise OMFITexception(f"{btip_nml} ip_direction does not match g-file")
        else:
            printw(f" > Running with bt {btd} and ip {ipd} >> Be sure this is correct")

as well as logic updating the gpec.in ip and bt directions when loading an equilibrium file.

@jmhanson

Thanks, Nik. I think I see the problem now. The ip_direction and bt_direction inputs are set correctly in coil.in. However, we have coil_flag = .false. in gpec.in. So I think I now understand that for intrinsic EF calculations using SURFMN, we need to put ip_direction and bt_direction in gpec.in because we are keeping the coil.in file disabled.

[logan-nc]

@ysm1887

It seems like the resonant coupling matrix looks the same.

@jmhanson

After looking at the documentation again I think this makes sense. And I do see a difference when I map to real space, as expected.I think we should also see that helicity makes a difference when we try to optimize I-coil corrections for the EF response.

The "working coordinate system" inside of DCON uses a sign convention such that the resonant surfaces are always positive m. So @ysm1887, the resonant coupling matrix is unchanged by helicity and the "dominant mode" that is most resonant will always peak in positive m.

@jmhanson, the EF from surfmn can get messed up if we do not get the helicity right. You can see it plays a role in how GPEC reads and converts the surfmn harmonics to the "working coordinate system" from DCON (it flips the order in which they are looped from the default mmin,mmax,1 to mmax,mmin,-1 where the 3rd number is the step direction.

GPEC/gpec/gpout.f

Lines 1126 to 1139 in c2b5ca4

	IF (data_flag) THEN
	IF(verbose) WRITE(*,*)"Reading external fields on the "//
	$ "control surface"
	ALLOCATE(dcosmn(mmin:mmax,nmin:nmax),
	$ dsinmn(mmin:mmax,nmin:nmax),rawmn(mmin:mmax,nmin:nmax))
	IF (data_type=="surfmn") THEN
	scale=1e-4
	IF (helicity<0) THEN
	i1=mmax
	i2=mmin
	i3=-1
	ENDIF
	ENDIF
	CALL ascii_open(in_unit,ifile,"old")

so I suppose that if the helicity is positive then we get lucky (the default helicity in GPEC is +1.0) and it all works out, but the EF will have the wrong helicity if the plasma current and bt are in opposite directions and we do not specify that explicitly in the namelists. Note @parkjk wrote all this interface before I was involved, so I don't have any plots to show proving this works out correctly for both cases. If you make some (always good to double check!) please, post them to this issue for future reference.