COCOS conversion not working, fix implemented

[henrikjaerleblad]

The conversion between different COCOS is not working in the current ascot5 version. This is important when loading eqdsk files with a COCOS that ASCOT cannot use, i.e. any COCOS other than 3.

I implemented a fix by adding a fix in a5py/templates/importdata.py (line 228) and a new COCOS class and functions in a5py/physlib/cocos.py

def fromCocosNtoCocosM(eqd, cocos_m, phiclockwise=None, weberperrad=None):
    """Transform equilibrium from cocos_n (determined from eqd, see below) to cocos_m.

    Parameters
    ----------
    eqd : dict
        Dictionary from reading the EQDSK file.
    cocosm : int
        Target COCOS.

    Returns
    -------
    eqdout : dict
        Equilibrium data converted to cocos_m.
    """

    cocos_n = assign(eqd["qpsi"][0], eqd["cpasma"], eqd["bcentr"], eqd["simagx"], eqd["sibdry"], phiclockwise, weberperrad)

    transform_dict = transform_cocos(cocos(cocos_n), cocos(cocos_m))

    # Define output
    eqdout = copy.deepcopy(eqd)
    # eqdout["nx"] = eqd["nx"] # For clarity
    # eqdout["ny"] = eqd["ny"] # -||-
    eqdout["rdim"] = eqd["rdim"]*transform_dict["R"]
    eqdout["zdim"] = eqd["zdim"]*transform_dict["Z"]
    eqdout["rcentr"] = eqd["rcentr"]*transform_dict["R"]
    eqdout["rleft"] = eqd["rleft"]*transform_dict["R"]
    eqdout["zmid"] = eqd["zmid"]*transform_dict["Z"]
    eqdout["rmagx"] = eqd["rmagx"]*transform_dict["R"]
    eqdout["zmagx"] = eqd["zmagx"]*transform_dict["Z"]
    eqdout["simagx"] = eqd["simagx"]*transform_dict["PSI"]
    eqdout["sibdry"] = eqd["sibdry"]*transform_dict["PSI"]
    eqdout["bcentr"] = eqd["bcentr"]*transform_dict["B"]
    eqdout["cpasma"] = eqd["cpasma"]*transform_dict["I"]
    eqdout["fpol"] = eqd["fpol"]*transform_dict["F"]
    eqdout["pres"] = eqd["pres"]*transform_dict["PRES"]
    eqdout["ffprime"] = eqd["ffprime"]*transform_dict["FFPRIME"]
    eqdout["pprime"] = eqd["pprime"]*transform_dict["PPRIME"]
    eqdout["psi"] = eqd["psi"]*transform_dict["PSI"]
    eqdout["qpsi"] = eqd["qpsi"]*transform_dict["Q"]
    #eqdout["psiaxis"] = psiaxis # Wrong key? psiaxis is called simagx according to freeqdsk docs
    #eqdout["psiedge"] = psiedge # Wrong key? psiedge is called sibdry according to freeqdsk docs
    #eqdout["nbdry"] = eqd["nbdry"]
    #eqdout["nlim"] = eqd["nlim"]
    eqdout["rbdry"] = eqd["rbdry"]*transform_dict["R"]
    eqdout["zbdry"] = eqd["zbdry"]*transform_dict["Z"]
    eqdout["rlim"] = eqd["rlim"]*transform_dict["R"]
    eqdout["zlim"] = eqd["zlim"]*transform_dict["Z"]

    return eqdout

def transform_cocos(cc_in: COCOS, cc_out: COCOS, sigma_Ip = None, sigma_B0 = None, ld = (1, 1), lB = (1, 1), exp_mu0 = (0, 0)):
    """
    Returns a dictionary of the multiplicative factors to transform COCOS from `cc_in` to `cc_out`

    Parameters
    ----------
    sigma_Ip : Union[Tuple[int, int], None]
        A tuple of the (Input, Output) current sign or nothing
    sigma_B0 : Union{NTuple{2,Int},Nothing}` - A tuple of the (Input, Output) toroidal field sign or nothing
    ld : NTuple{2,<:Real}
        A tuple of the (Input, Output) length scale factor. Default = (1,1)
    lB : NTuple{2,<:Real}
        A tuple of the (Input, Output) magnetic field scale factor. Default = (1,1)
    exp_mu0 : NTuple{2,<:Real}
        A tuple of the (Input, Output) mu0 exponent (0, 1). Default = (0,0)

    Returns
    -------
    transforms : dict
        Transform multiplicative factors to be able to convert from `cc_in to `cc_out`
    """

    ld_eff = ld[1] / ld[0]
    lB_eff = lB[1] / lB[0]
    exp_mu0_eff = exp_mu0[1] - exp_mu0[0]

    sigma_RpZ_eff = cc_in.sigma_RpZ * cc_out.sigma_RpZ

    if sigma_Ip is None:
        sigma_Ip_eff = cc_in.sigma_RpZ * cc_out.sigma_RpZ
    else:
        sigma_Ip_eff = sigma_Ip[0] * sigma_Ip[1]

    if sigma_B0 is None:
        sigma_B0_eff = cc_in.sigma_RpZ * cc_out.sigma_RpZ
    else:
        sigma_B0_eff = sigma_B0[0] * sigma_B0[1]

    sigma_Bp_eff = cc_in.sigma_Bp * cc_out.sigma_Bp
    exp_Bp_eff = cc_out.exp_Bp - cc_in.exp_Bp
    sigma_rhotp_eff = cc_in.sigma_rhotp * cc_out.sigma_rhotp

    mu0 = 4 * 3.14159265358979323846 * 1e-7  # pi is used directly for more precision

    transforms = {}
    transforms["R"] = ld_eff
    transforms["Z"] = ld_eff
    transforms["PRES"] = (lB_eff ** 2) / (mu0 ** exp_mu0_eff)
    transforms["PSI"] = lB_eff * (ld_eff ** 2) * sigma_Ip_eff * sigma_Bp_eff * ((2 * 3.14159265358979323846) ** exp_Bp_eff) * (ld_eff ** 2) * lB_eff
    transforms["TOR"] = lB_eff * (ld_eff ** 2) * sigma_B0_eff
    transforms["PPRIME"] = (lB_eff / ((ld_eff ** 2) * (mu0 ** exp_mu0_eff))) * sigma_Ip_eff * sigma_Bp_eff / ((2 * 3.14159265358979323846) ** exp_Bp_eff)
    transforms["FFPRIME"] = lB_eff * sigma_Ip_eff * sigma_Bp_eff / ((2 * 3.14159265358979323846) ** exp_Bp_eff)
    transforms["B"] = lB_eff * sigma_B0_eff
    transforms["F"] = sigma_B0_eff * ld_eff * lB_eff
    transforms["I"] = sigma_Ip_eff * ld_eff * lB_eff / (mu0 ** exp_mu0_eff)
    transforms["J"] = sigma_Ip_eff * lB_eff / ((mu0 ** exp_mu0_eff) * ld_eff)
    transforms["Q"] = sigma_Ip_eff * sigma_B0_eff * sigma_rhotp_eff

    return transforms

The fix is based on the equations that describe how to transform between any COCOS in O. Sauter et al, CPC, 2013. I will create a new branch locally and submit a pull request to fix this issue, in line with the developer docs found in https://ascot4fusion.github.io/ascot5/developing.html.

Please let me know if you have any comments!

[henrikjaerleblad]

I should add that the current implementation of the fix can convert from any COCOS to odd COCOS. I am currently investigating why. Nevertheless, ASCOT needs COCOS 3, so the current fix can still be used to convert any eqdsk file with any COCOS to the COCOS required by ASCOT

[henrikjaerleblad]

I found it! It was just a matter of putting the phiclockwise keyword argument in cocos.py/assign() to True instead of False.

In short, the current fix thus works to convert from any COCOS to any COCOS :)

(I verified via the cocos.py/assign() function)