Skip to content

ExB/electron profile conversion subtracts the diamagnetic rotation #104

Description

@krystophny

Summary

calc_omega_profile documents that iread_omega_ExB adds the ion diamagnetic term to obtain ion rotation, but the implementation subtracts it. This conflicts with maintained M3D-C1 post-processing definitions, the published M3D-C1 rotation relation, and a native initialization-only discriminator.

This report concerns the relative sign inside M3D-C1. It does not infer a foreign equilibrium's COCOS and does not propose flipping any plotted result.

Current source

At public physics revision c7f9c14a26fc72dc679598609b2c424952b0300c, unstructured/gradshafranov.f90 has:

else if(iread_omega_ExB.ne.0) then
   ! we're reading in ExB rotation; add ion diamagnetic term
   ...
   dia = db*(ppval/nval - z_ion*(1.+thermal_force_coeff)*tep &
        - z_ion*te0*np/nval)
...
omega_spline%y(i) = omega_spline_0%y(i) - dpsii*dia

The iread_omega_e path uses the same final subtraction.

Independent M3D-C1 contracts

The bundled diagnostics consistently use the opposite relation:

  • unstructured/idl/plot_omega.pro: w_star_i = db*piprime/den, then omega_ExB = v_omega - w_star_i;
  • unstructured/idl/read_field.pro: omega_*i = db*d(p-p_e)/dpsi/den, then omega_ExB = omega - w_star_i;
  • unstructured/idl/rmp.pro and schaffer_plot.pro use omega_ExB = omega_i - w_star_i.

These all imply omega_i = omega_ExB + omega_*i.

Ferraro et al., Nuclear Fusion 53 (2013) 073042, DOI 10.1088/0029-5515/53/7/073042, printed page 5, states for M3D-C1 that the perpendicular rotations satisfy omega_{e,i}=omega_ExB+omega_{*e,*i}, with omega_*i=(d p_i/dpsi)/(n_e e).

In M3D-C1's explicit right-handed cylindrical representation,
B=grad(psi) x grad(phi)+F grad(phi), E=-Phi'(psi) grad(psi), and positive-ion v_*i=B x grad(p_i)/(q_i n_i B^2). The toroidal coefficients again give omega_i=omega_ExB+(1/(q_i n_i)) dp_i/dpsi.

Native discriminator

An initialization-only run selected iread_omega_ExB=1, db=0.0322, and wrote the native profile before any time step. GDB recorded dpsii=3.1702203540215126 after define_profiles. At normalized psi zero:

normalized input omega_ExB       0.06959516098919606
native dia                       -0.013683542009059833
native omega_i                   0.1129750044
current minus prediction         0.11297500438142596
documented plus prediction       0.02621531759696616
native/current residual          1.86e-11
native/documented separation     8.68e-2

Thus the executable exactly follows the final minus sign; this is not a display, interpolation, or COCOS ambiguity.

Suggested resolution

Please confirm the intended iread_omega_ExB and iread_omega_e contracts. If the comments, bundled diagnostics, and published relation are intended, change the final operation to addition and add a regression with nonzero pressure gradients. The regression should cover both input paths, both signs of the diamagnetic term, and the intended behavior of iflip_j/iflip_v rather than assuming those modifiers.

No M3D-C1 production source was changed during this audit. The reproducer, hashes, source-token checks, raw GDB trace, and exact arithmetic are maintained externally in verify_m3dc1_profile_rotation.py.

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