diff --git a/CLAUDE.md b/CLAUDE.md index 9058abeae..428772c67 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -257,7 +257,7 @@ decomposition exp(im*theta - in*phi). A comprehensive reference is in - **helicity** = ipd * btd (+1 = RH, -1 = LH), computed in `gpec_main` (`gpec/gpec.f`) - **nn** (toroidal mode number) is always positive; resonant **m** is always positive - **F = R*Bt** is forced positive via ABS() in `read_eq_efit` (`equil/read_eq.f`) -- **q** is NOT forced positive (read directly from EFIT) +- **q** is recomputed by field-line integration (`direct_run` in `equil/direct.f`; `inverse_run` in `equil/inverse.f` for inverse formats); the input file's q profile is unused, and the recomputed q is always positive - **omega_E** is positive for rotation in the direction of the toroidal coordinate zeta - The code does NOT use the COCOS standard - For SURFMN interface: `m_surfmn = helicity * m_gpec` diff --git a/docs/sign_conventions.rst b/docs/sign_conventions.rst index 42227a0b1..3754c615e 100644 --- a/docs/sign_conventions.rst +++ b/docs/sign_conventions.rst @@ -96,9 +96,9 @@ Helicity is computed in the ``gpec_main`` program (``gpec/gpec.f``): - ``ip_direction`` and ``bt_direction`` are set in ``coil.in``. "positive" means CCW viewed from above; "negative" means CW. -- **helicity = +1**: right-handed (RH) --- :math:`B_t` and :math:`I_p` in the +- **helicity = +1**: right-handed (RH), :math:`B_t` and :math:`I_p` in the same direction. -- **helicity = -1**: left-handed (LH) --- :math:`B_t` and :math:`I_p` opposed. +- **helicity = -1**: left-handed (LH), :math:`B_t` and :math:`I_p` opposed. - The helicity value is stored in the ``gpec_control_output`` netcdf file. @@ -121,12 +121,22 @@ Safety Factor :math:`q` ======================== - Defined as :math:`q = (\mathbf{B} \cdot \nabla\zeta) / (\mathbf{B} \cdot \nabla\theta)`. -- **Not forced positive.** Read directly from the EFIT g-file without sign - manipulation (unlike :math:`F`). -- In standard tokamak operation, EFIT provides :math:`q > 0`. -- Can be adjusted via ``newq0`` in ``equil.in``, which modifies :math:`F` - to match while preserving the Grad-Shafranov solution - (``direct_run`` in ``equil/direct.f``). +- **Recomputed, not read.** For direct equilibria (EFIT g-files), the q + profile in the file is read but never used. ``direct_run`` in + ``equil/direct.f`` rebuilds q on each flux surface from a field-line + integration. +- **Always positive for g-file input.** The integration uses the + sign-normalized :math:`\psi` map and :math:`|F|`, so q comes out positive + for any combination of signs in the g-file. Reversing all signed + quantities in a g-file (:math:`\psi`, ``simag``, ``sibry``, ``cpasma``, + q) reproduces the original DCON results bit for bit. +- **Inverse equilibria recompute q too.** The main inverse path + (``inverse_run`` in ``equil/inverse.f``, used by the CHEASE, JSOLVER, + TRANSP, and analytic formats) rebuilds q on each surface by a flux-surface + integration, exactly as the direct path does; the input file's q profile + is read into ``sq_in`` but not used, and the recomputed q is positive. The + one exception is the ``chease4`` format (``inverse_chease4_run``), which + takes q directly from the input file. Mode Numbers :math:`m` and :math:`n` @@ -157,7 +167,8 @@ Why Positive :math:`m` Is Always Resonant Resonant surfaces are found by ``sing_find`` in ``dcon/sing.f``. It performs a binary search for flux surfaces where :math:`m = n \cdot q`. -Since :math:`n > 0` (by convention) and :math:`q > 0` (standard tokamak): +Since :math:`n > 0` (by convention) and :math:`q > 0` (the recomputed q is +always positive, see the safety factor section above): .. math:: @@ -284,6 +295,19 @@ system. The conventions described in this document are GPEC-native and predate COCOS. Users interfacing with COCOS-aware codes must manually translate between conventions. +For g-file input the sign part of the COCOS choice cannot matter: +``read_eq_efit`` sign-normalizes the :math:`\psi` map, forces +:math:`F > 0`, and discards the file's q profile, so g-files that differ +only in sign conventions produce identical results. Field and current +directions enter solely through ``ip_direction`` and ``bt_direction`` in +``coil.in``. The main inverse path likewise recomputes q rather than +inheriting it from the input file. What does matter is the flux +normalization: GPEC assumes the +g-eqdsk standard of poloidal flux per radian (Wb/rad). A file carrying the +full flux (the COCOS 11-18 family) yields a wrong q magnitude. Reports +that g-files must be supplied in one specific COCOS trace to this unit +requirement, since the sign choices are normalized away on read. + Quick Reference =============== @@ -308,8 +332,8 @@ Quick Reference - Always positive - Yes, ABS() * - :math:`q` (safety factor) - - From EFIT, typically positive - - No + - Recomputed by field-line integration + - Always positive * - :math:`n` (toroidal mode) - Always positive - By convention @@ -330,6 +354,8 @@ Source Code References - **Poloidal flux sign**: ``read_eq_efit`` in ``equil/read_eq.f`` - **F = R*Bt forced positive**: ``read_eq_efit`` in ``equil/read_eq.f`` - **q definition**: ``dcon/README`` +- **q from field-line integration**: ``direct_run`` in ``equil/direct.f`` +- **q recomputed (inverse path)**: ``inverse_run`` in ``equil/inverse.f`` - **Helicity computation**: ``gpec_main`` program in ``gpec/gpec.f`` - **ip/bt direction**: ``input/coil.in`` - **Poloidal mode range**: ``dcon`` program in ``dcon/dcon.F`` diff --git a/equil/direct.f b/equil/direct.f index 5842919cb..e97ab45e6 100644 --- a/equil/direct.f +++ b/equil/direct.f @@ -195,7 +195,6 @@ SUBROUTINE direct_run sq%name=" sq " sq%title=(/" psi "," f "," p "," q "/) q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) - IF(newq0 == -1)newq0=-q0 c----------------------------------------------------------------------- c revise q profile. c----------------------------------------------------------------------- @@ -204,7 +203,7 @@ SUBROUTINE direct_run f0fac=f0**2*((newq0/q0)**2-1) q0=newq0 DO ipsi=0,mpsi - ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0) + ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac diff --git a/equil/equil.f b/equil/equil.f index f471f404f..00ad1aad0 100644 --- a/equil/equil.f +++ b/equil/equil.f @@ -55,6 +55,11 @@ SUBROUTINE equil_read(unit, op_psihigh, op_psilow) READ(UNIT=in_unit,NML=equil_output) ENDIF CALL ascii_close(in_unit) +c----------------------------------------------------------------------- +c guard against negative newq0 (negative q is not supported). +c----------------------------------------------------------------------- + IF(newq0 < 0)CALL program_stop("newq0 < 0 is not supported: "// + $ "q must be positive. Use newq0 = 0 (default) or newq0 > 0.") IF(PRESENT(op_psihigh))THEN psihigh = op_psihigh IF(verbose) WRITE(*,*) "Reforming equilibrium with new psihigh" diff --git a/equil/inverse.f b/equil/inverse.f index d41b379ac..23405ea38 100644 --- a/equil/inverse.f +++ b/equil/inverse.f @@ -237,7 +237,6 @@ SUBROUTINE inverse_run ENDDO CALL spline_fit(sq,"extrap") q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) - IF(newq0 == -1)newq0=-q0 c----------------------------------------------------------------------- c revise q profile. c----------------------------------------------------------------------- @@ -246,7 +245,7 @@ SUBROUTINE inverse_run f0fac=f0**2*((newq0/q0)**2-one) q0=newq0 DO ipsi=0,mpsi - ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0) + ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac @@ -579,7 +578,11 @@ SUBROUTINE inverse_chease4_run(jaceq) ENDDO CALL spline_fit(sq,"extrap") q0=sq%fs(0,4)-sq%fs1(0,4)*sq%xs(0) - IF(newq0 == -1)newq0=-q0 +c----------------------------------------------------------------------- +c enforce positive q (chease4 takes q directly from the input file). +c----------------------------------------------------------------------- + IF(MINVAL(sq%fs(:,4)) <= 0)CALL program_stop + $ ("chease4 input q must be positive (GPEC requires q > 0).") c----------------------------------------------------------------------- c revise q profile. c----------------------------------------------------------------------- @@ -588,7 +591,7 @@ SUBROUTINE inverse_chease4_run(jaceq) f0fac=f0**2*((newq0/q0)**2-one) q0=newq0 DO ipsi=0,mpsi - ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0) + ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac diff --git a/equil/read_eq.f b/equil/read_eq.f index b5e568a79..76793c06d 100644 --- a/equil/read_eq.f +++ b/equil/read_eq.f @@ -1498,13 +1498,12 @@ SUBROUTINE read_eq_dump c----------------------------------------------------------------------- c revise q profile. c----------------------------------------------------------------------- - IF(newq0 == -1)newq0=-q0 IF(newq0 /= 0)THEN f0=sq%fs(0,1)-sq%fs1(0,1)*sq%xs(0) f0fac=f0**2*((newq0/q0)**2-one) q0=newq0 DO ipsi=0,mpsi - ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2)*SIGN(one,newq0) + ffac=SQRT(1+f0fac/sq%fs(ipsi,1)**2) sq%fs(ipsi,1)=sq%fs(ipsi,1)*ffac sq%fs(ipsi,4)=sq%fs(ipsi,4)*ffac rzphi%fs(ipsi,:,3)=rzphi%fs(ipsi,:,3)*ffac