Expose per-point qflux current-channel contributions - #151
Conversation
isw_qflux_NA=1 with in_file_pert='none' assembles an identically zero non-axisymmetric source. That zero vector is the Arnoldi seed, so its normalization divides by zero (0/0 -> NaN) and zhseqr then fails to converge at every step, spinning without qflux output. Reject the configuration in check() with a clear message instead.
Add an optional right-hand side to the non-axisymmetric solve: the drive of a single-helicity perturbation given by a radial corrugation field (delta B^s/B_0^phi)_mn and an electrostatic potential harmonic Phi_mn. The vpar corrugation piece enters the A1/A2 force channels with velocity weights a3m and x^2 (sigma-odd), the ExB piece with weights x^-1 and a1m (sigma-even); the two new velocity moments are computed in compute_sources. Input via new &ntv_input variables (isw_hel_drive, m_theta_hel, hel_brad_re/im, hel_phim_re/im), default off. Band profiles, moments and dimensional prefactors follow a CAS-checked derivation (see issue #118). Off by default: golden-record QL summary.h5 is reproduced to 1e-15 against a main build.
Review findings on the drive commit: fold the x^(3k-4) weight into a single power so the x=0 quadrature node stays finite; skip the moment quadrature on the precomputed-matrix-element path and reject isw_hel_drive at run time when the moments are unavailable; reject the unvalidated relativistic branch, reconstruction mode and unset species data; deallocate asource_hel in collop_unload.
The zero-source guard added upstream rejects isw_qflux_NA=1 with in_file_pert='none'. With the misalignment drive, that configuration is valid: the drive supplies the non-axisymmetric source. Relax the guard so it fires only when the drive is off or all its amplitudes are zero.
The opposite-helicity drive (m,n) -> (m,-n) enters add_helical_source as the complex conjugate of the band profile (phase convention exp(+i(m theta + n phi))). The assembly stencil is a real-coefficient linear map, so it must conjugate the assembled source pointwise; a sign or phase error in the stencil breaks this. The existing checks use a profile whose real and imaginary parts vary independently, so they do not exercise conjugate symmetry. Add a check that conjg(profile) yields conjg(source) to 1e-14, plus a vacuity guard that the source carries an imaginary part. Verification: - ctest helical_source_test: Passed (0.07 s) - ctest -LE external-data: 18/18 passed - red-before: injecting a complex coefficient (1,0.1) into one stencil term makes the new check error stop with 'conjugate profile did not conjugate the source pointwise' (exit 1)
Decompose the phi-integrated current channel qflux(2,k) = sum_col flux_vector(2,col)*source_vector(col,k) into its per-field-line-point contributions in the mode-1 local (axi) solver, which HELCORE uses. The new pointwise_current_mod reuses the already-verified flux_vector current row: its counter-passing block enters with the explicit minus sign (parallel current is odd under sigma -> -sigma) and the step_factor phi-quadrature weight is baked in, so the decomposition introduces no new sign, weight, or normalization. flux_vector's columns partition into disjoint per-istep blocks, so the block dot products sum back to qflux(2,k) exactly. The solver writes current_profile.dat (phi_mfl, per-drive contributions, and current densities) guarded by isw_hel_drive, so runs without the helical drive are bitwise unchanged. A hard self-check stops the run if the field-line sum does not reproduce qflux(2,k). Verification: - 17/17 TEST ctests pass; new pointwise_current_test added. - Red-before: breaking the block partition (nblk without the npassing+1) gives "ERROR STOP FAIL: per-point contributions do not sum to the qflux current channel"; restored -> Passed. - Local HELCORE simple_l1 run, all three drive cases (nminus_phi, nminus_phi_im, nplus_phi): gamma_current/D31/D32 reproduce the accepted aCluster response (-7.338706/-20.944508, -2.711060/-3.138948, -2.369674/-1.825414) to ~1e-11; each current_profile.dat reconstructs qflux(2,k) to ~1e-15 relative over 961 finite points.
krystophny
changed the title
![slopqueue[bot]](https://avatars.githubusercontent.com/u/48486978?s=60&v=4){kind=small}
There was a problem hiding this comment.
Review verdict: Request changes
Summary
This PR adds a single-helicity misalignment drive ("helical drive") to the NEO-2-QL non-axisymmetric solver: new physics modules (helical_source_mod, helical_response_mod, qflux_profile_mod, join_diagnostics_mod), collision-operator endpoint fixes, NTV output geometry synchronization, and extensive new tests. The physics and I/O structure are sound, but one production-code boundary bug needs fixing before merge.
Findings
-
major
NEO-2-QL/helical_source_mod.f90:70-80— Inadd_forward_source, the loopdo i = ibeg, iendenters themod(i-ibeg,2)==1branch wheni = iend(assumingmod(iend-ibeg,2)==1), which accessesnpl(i+1)=npl(iend+1). Sincenplis dimensionedibeg:iend, this is an out-of-bounds read (undefined behavior — crash or silent garbage). The siblingadd_backward_sourcecorrectly guards its right boundary withif (i >= iend) cycle, butadd_forward_sourceonly guards the left (if (i <= ibeg) cycle) and never guards the right. Fix: addif (i >= iend) cycle(or equivalent boundary guard) inadd_forward_sourcebefore themod(i-ibeg,2)branch, mirroring the backward source's boundary protection. -
minor
NEO-2-QL/ntv_mod.f90:240-246—has_perturbation_filechecks for'none','None','NONE'as separate comparisons rather than a single case-insensitive check. A user entering'NONE '(with trailing space) afterADJUSTLwould still match viaTRIM, but any other casing (e.g.,'nOne') would not. Considerto_loweror a singlescanfor robustness. Non-blocking. -
minor
NEO-2-QL/propagator.f90:~965— The qflux-contribution validation comparesABS(cp_channel - qflux(2,k))against1d-10 * MAX(1, ABS(qflux(2,k))). For large numbers of field-line points, accumulated floating-point round-off in the sum-of-products may exceed this tolerance, causing a falseERROR STOPon some platforms. Consider relaxing to1d-8or1d-6relative. Non-blocking.
Verdict
Request changes — the out-of-bounds access in add_forward_source (finding 1) is undefined behavior in production code and must be fixed before merge. The remaining findings are minor and can be addressed at the author's discretion.
There was a problem hiding this comment.
Review verdict: Request changes
Summary
PR #151 adds a single-helicity misalignment ("helical") drive to the NTV non-axisymmetric solver, alongside new modules (helical source/response, qflux profile, join diagnostics), collision-operator zero-endpoint fixes, phi-divider parity refinement, and expanded tests.
Findings
-
major
NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90:~3196-3260— The single-helicity drive block has unbalanced parentheses (e.g.hel_phase_fac = imun*(DBLE(m_theta_hel)*aiota + DBLE(m_phi))has an extra closing paren, andrho_hel = .../(bmod0*1.0d4))has a spurious)). This will not compile. Required fix: correct all parenthesization in the helical drive section. -
major
NEO-2-QL/ripple_solver_ArnoldiOrder2_test.f90:~4261— Thesource_fluxsubroutine callsadd_helical_sourcefour times, but the first two and second two are duplicate pairs (same argumentsq_hel_b/asource(0:lag,2)andq_hel_b/asource_hel(0:lag,2), etc.). This double-counts the helical drive, corrupting the physics. Required fix: remove the duplicate calls or verify whether both A1/A2 channels are intended. -
minor
TEST/test_helical_response.f90:~24—if (any(gamma_current /= [2.0_dp, -3.0_dp]))has an extra)(two closing parens for oneany(open). This will not compile. Fix:any(gamma(3,1:2) /= [2.0_dp, -3.0_dp]).
Verdict: Request changes — the unbalanced parentheses will break compilation, and the duplicated source calls produce incorrect physics. Fix these before merge.
1 participant
Exposes a per-field-line-point partition of the raw mode-1 local-solver
qflux(2,k)current channel.Stack: top of the open flux-pumping response stack. The parent is
feat/flux-pumping/source-conjugacyat86f9503. No merge is requested duringthe integration freeze.
The diagnostic partitions
flux_vectorinto its existing per-step columnblocks. Their dot products sum to
qflux(2,k)exactly. The output isqflux_current_contributions.dat; it contains raw normalized quadraturecontributions, not pointwise physical parallel-current density.
The correction at
fd22f71removes the former division by centered toroidalangle spacing. That spacing was not the solver's direction-dependent
Runge-Kutta/Simpson measure and could not produce a local current. Physical
current remains a later normalization and force-application step.
Invariant: the kinetic operator, source, discretization, convergence criteria,
current-row sign, normalization, units, coordinate conventions, boundary
conditions, ABI, and memory layout are unchanged. Runs without the helical
drive remain bitwise unchanged. A hard guard still rejects any partition that
does not sum to
qflux(2,k).Verification
Test fails on main
mainhas no qflux-profile diagnostic. The rejected behavior is present onthe previous PR head:
This divided a weighted qflux contribution by an unrelated centered angle
spacing and labelled the result as physical current density.
Test passes after fix