KineticForces - FEATURE - Rotation prediction: a momentum-balance module fed by tabulated torque surfaces #461

Why

Several questions we want the code to answer are the same computation seen from different sides:

• What steady-state rotation do we predict if these 3D fields are on, and what if they are switched off? Today the answer stops at "the NTV torque at the given rotation profile".
• What is the true neoclassical offset rotation profile for a given perturbation, rather than the rough ~2·ω_*T estimate?
• How much error field can a correction coil cancel before its own non-resonant torque costs the rotation that keeps the penetration threshold up (the NTV-limits analysis in ErrorFields; its current form assumes torque ∝ rotation far from the offset and a positive torque).
• Profile-consistent critical resonant field (Tearing.CriticalResonantField, Tearing - FEATURE - Add critical resonant field calculations to the tearing workflow #415, SLAYER based b_crit calculations need to be ported #371, Tearing.CriticalResonantField - Improve pole avoidance, input options, and grid refinement #416): the layer's rotation should come from a momentum balance that includes the 3D-field torques, not only from the unperturbed profile.
• Soon: the resonant J×B torque from the resistive perturbed equilibrium as a second torque source (PE - BUG - Im(py) is ~900x larger than Fortran GPEC, producing a spurious toroidal torque #401, PE toroidal torque changes sign under grid refinement while its magnitude stays stable to ~1% #370, PE toroidal torque nonzero for ideal-MHD runs — track imag/real ratio as a Hermiticity check #333 track its sign and Hermiticity problems).

All of these are a momentum balance with torque sources that depend on the rotation they act on. The torque sources are expensive (a kinetic torque evaluation is minutes); the balances are cheap. If a transport or balance solver calls the kinetic solve at every iteration or time step the whole thing is unusable, so the framework has to separate the two from the start.

Principle: torque surfaces, evaluated once, consumed by any solver

Every torque source publishes the same object:

TorqueSurface
  psi          :: Vector       # kinetic ψ_N grid the profile is resampled onto
  rotation     :: Vector       # rotation-shift grid Δω (rad/s), adaptive, contains 0
  torque       :: Matrix       # cumulative torque T(ψ, Δω) [N·m], per unit perturbation amplitude²
  dtorque_dpsi :: Matrix       # torque density
  source, method, perturbation label, amplitude normalization, sign convention

built once per perturbation by scanning the rotation with an adaptive sampler (the NTV torque has sharp features: the superbanana-plateau resonance near ω_E ≈ 0, bounce-harmonic resonances, the offset crossing; a fixed grid misses them or wastes evaluations). Solvers interpolate the surface and never call the physics kernel. Sources add: NTV now, resonant EM later, each its own surface with the same contract. Amplitude scaling is quadratic, so a surface per unit current serves any current.

The kinetic torque already returns the ψ-resolved cumulative torque (torque_profile) and the KineticProfileSplines container has a vector constructor, so a rigid shift of ω_E with the diamagnetic terms held fixed is a few lines; that is the scan variable for NTV (braking changes E×B, not ω_*). The sign convention between T_φ and the rotation it acts on must be a named constant with a runtime check (a restoring torque has negative slope at its zero crossing), never an assumption.

Proposed home

A new module after KineticForces and Tearing and before Analysis, working name MomentumBalance (or Transport; naming is open). It owns the balance solvers and the TorqueSurface type; it does not compute torques. Pieces that belong elsewhere:

• Utilities.AdaptiveSampling: the generic adaptive 1-D sampler (vector-valued function, initial grid, refinement on curvature and sign change, deterministic, point cap). Reusable for every rotation scan and for any other sharp-feature scan.
• Equilibrium.KineticProfiles: shift_exb_rotation(profiles, Δ).
• Momentum inputs: the kinetic HDF5 schema already carries omega_tor, chi_e, chi_phi, and Tearing - FEATURE - Add critical resonant field calculations to the tearing workflow #415 added viscous_input_type / viscous_input (profile, scalar, or Prandtl number) to the Tearing control, with reviewers asking for a scalar τ_φ interface. Propose one resolver shared by Tearing.CriticalResonantField and MomentumBalance: profile column → scalar χ_φ → scalar τ_φ (converted to a flat χ_φ from the momentum content) → fallback to χ_e, so both consumers see the same χ_φ(ψ).
• Outputs under a MomentumBalance/ group in gpec.h5 (predicted rotation profiles, offset profile, the torque surfaces with their provenance), annotated per the HDF5 conventions.

Initial capabilities, ranked

1. 0-D balance from a torque surface (in progress on the ErrorFields NTV-limits branch): Δ·T_0/ω_ref = T(Δ)·I² solved for the rotation shift; threshold scaling with an exponent; loss of the root as the braking bifurcation. Acceptance: reduces exactly to the linear-budget model when the torque is constant; the DIII-D error-field example shows the torque crossing zero inside the scanned span.
2. Neoclassical offset profile: the zero of T′(ψ, Δ) in Δ at each ψ, with the rough ~2·ω_*T estimate reported next to it. Acceptance: the offset lies inside the scan span everywhere the torque density is non-negligible.
3. Steady 1-D momentum balance: −(1/V′) d/dψ [V′ n m ⟨R²⟩ χ_φ |∇ψ|² dω/dψ] = Σ sources with the NTV surface as source, χ_φ from the shared resolver, sources from an input torque profile or the unperturbed rotation, and the edge rotation as boundary condition. Outputs the predicted ω(ψ) with and without the 3D field. Acceptance: with zero 3D field it returns the input profile; with a surface that is linear in Δ it matches the analytic solution.
4. Profile-consistent b_crit: feed the balanced rotation at each rational surface to CriticalResonantField.
5. Resonant EM torque surface from the resistive PE, same contract, added as a second source; needs PE - BUG - Im(py) is ~900x larger than Fortran GPEC, producing a spurious toroidal torque #401/PE toroidal torque changes sign under grid refinement while its magnitude stays stable to ~1% #370/PE toroidal torque nonzero for ideal-MHD runs — track imag/real ratio as a Hermiticity check #333 settled first.
6. Time-dependent balance (ramps, bifurcation dynamics) once 3 exists; still interpolating the surfaces, never re-evaluating the kernel per step.

Open questions

• Which rotation is the balance variable: rigid shift of ω_E with ω_* fixed (proposed), or a scale factor, or ω_φ with the diamagnetic terms re-derived from the balanced profiles?
• ω_ref for the 0-D friction term: momentum-weighted mean of ω_φ (proposed) versus the rotation at the dominant rational surface.
• Boundary condition and edge/SOL treatment for the 1-D solve (Physics based NTV edge bound is needed #315 is the NTV edge bound).
• The resonant torque depends on rotation through the layer response, not through a rigid shift; its scan variable may differ.
• Should the surfaces be per unit current (ErrorFields use) or per unit boundary field amplitude (general use)? Both are one scaling away.

Related: #415, #371, #416, #423 (ψ-resolved torque response profiles), #315, #401, #370, #333.