Skip to content

Sign reference disagrees with PENTRC passing harmonic and co-current map #276

Description

@krystophny

The current sign-conventions reference at aaf720d475df6616935a56ea5914a68df10575d0 disagrees with executable source in two sign-sensitive places.

  1. It documents the passing-particle effective harmonic as ell_eff = ell - sigma*n*q:

Energy Integral Resonance
-------------------------
In the PENTRC energy integral (``xintgrnd`` in ``pentrc/energy.f90``), the
resonance denominator involves:
.. math::
n \omega_E + \ell_{\mathrm{eff}} \omega_b \sqrt{x} + n \omega_D x
where :math:`\omega_b` is the bounce frequency divided by :math:`x`,
:math:`\omega_D` is the magnetic precession frequency,
:math:`\ell_{\mathrm{eff}} = \ell - \sigma n q` is the effective bounce harmonic,
and :math:`x = E/T` is the normalized energy. The sign of :math:`\omega_E`
determines the direction of resonance in velocity space.

All three executable paths use the plus sign:

  • pitch.f90:

    GPEC/pentrc/pitch.f90

    Lines 239 to 247 in aaf720d

    if(lmda>bobmax)then
    nueff = nuk/(2*epsr)
    lnq = real(ell,r8)
    else
    nueff = nuk
    lnq = ell+n*q
    endif
    ! note: currently ignores -wb case for trapped
    xint = xintgrl_lsode(wn,wt,we,wd,wb,nueff,ell,lnq,n,psi,lmda,method,record_this)
  • pitch.f90 runtime integrand:

    GPEC/pentrc/pitch.f90

    Lines 357 to 369 in aaf720d

    ! energy integration of resonance operator
    if(x<=pitch_bobmax)then ! circulating particles
    nueff = pitch_nuk
    lnq = pitch_ell+pitch_n*pitch_q
    xint = xintgrl_lsode(pitch_wn,pitch_wt,pitch_we,wd,wb,nueff,pitch_ell,lnq,pitch_n,&
    pitch_psi,real(x,r8),pitch_method,op_record=.false.)
    xint = xint + xintgrl_lsode(pitch_wn,pitch_wt,pitch_we,wd,-wb,nueff,pitch_ell,lnq,pitch_n,&
    pitch_psi,real(x,r8),pitch_method,op_record=.false.)
    else ! trapped particles
    nueff = pitch_nuk/(2*pitch_epsr) ! effective collisionality
    lnq = real(pitch_ell,r8)
    xint = xintgrl_lsode(pitch_wn,pitch_wt,pitch_we,wd,wb,nueff,pitch_ell,lnq,pitch_n,&
    pitch_psi,real(x,r8),pitch_method,op_record=.false.)
  • torque.F90:

    GPEC/pentrc/torque.F90

    Lines 591 to 597 in aaf720d

    if(lmda>(bo/bmax)) then
    sigma = 0 !trapped
    else
    sigma = 1 !passing
    endif
    lnq = l+sigma*n*q

So the maintained source implements ell_eff = ell + sigma*n*q.

  1. The page says positive native omega_E is co-current for right-handed plasmas and counter-current for left-handed plasmas:

Rotation Velocity Conventions (PENTRC)
======================================
:math:`\omega_E` (E x B Rotation)
----------------------------------
- Column 6 of the PENTRC kinetic profile file: :math:`\omega_E` in rad/s.
- Read by the ``read_kin`` subroutine in ``pentrc/inputs.f90``.
- **Sign convention**: positive :math:`\omega_E` means rotation in the
direction of the toroidal coordinate :math:`\zeta`.
- Since :math:`\phi` direction depends on helicity, positive :math:`\omega_E`
is effectively **co-current for RH plasmas** and **counter-current for
LH plasmas**.

But the physical toroidal map in the coil field path is

phi_CCW = -helicity * (2*pi*zeta + delta_phi)

GPEC/coil/field.F

Lines 166 to 193 in aaf720d

! Loop over every theta, zeta point on the flux surface
DO itheta=1,cmtheta
CALL bicube_eval_external(crzphi,psi,ctheta(itheta),1,
$ spline_ipsi,spline_itheta,
$ crzphi_f,crzphi_fx,crzphi_fy)
rfac=SQRT(crzphi_f(1))
eta=twopi*(ctheta(itheta)+crzphi_f(2))
rr=cro+rfac*COS(eta)
zz=czo+rfac*SIN(eta)
jac=crzphi_f(4)
w11=(1+crzphi_fy(2))*twopi**2*rfac*rr/jac
w12=-crzphi_fy(1)*pi*rr/(rfac*jac)
delpsi=SQRT(w11**2+w12**2)
IF(wegt/=0)THEN
aspl_f_arr(itheta)=jac*delpsi
ELSE
aspl_f_arr(itheta)=1.0
ENDIF
w11=w11/delpsi
w12=w12/delpsi
DO izeta=1,cmzeta
phi=-helicity*(twopi*czeta(izeta)+crzphi_f(3))
xobs(itheta,izeta)=rr*COS(phi)
yobs(itheta,izeta)=rr*SIN(phi)
zobs(itheta,izeta)=zz

with helicity = ipd*btd and positive direction strings defined as CCW from above. Therefore the native +zeta sign relative to plasma current is

(-helicity) / ipd = -btd

(equivalently multiplication by ipd, since signs square to one). Co/counter-current depends on the absolute bt_direction, not helicity alone. For (ipd,btd,helicity)=(+1,+1,+1), positive native omega_E is clockwise and counter-current, contrary to the page.

These are documentation corrections only; this issue does not assert that the executable resonance or coordinate transforms should change.

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions