Skip to content

Expose reconstructed local current-row components#152

Open
krystophny wants to merge 27 commits into
mainfrom
feat/flux-pumping/local-parallel-response
Open

Expose reconstructed local current-row components#152
krystophny wants to merge 27 commits into
mainfrom
feat/flux-pumping/local-parallel-response

Conversation

@krystophny

Copy link
Copy Markdown
Member

Parent stack head: PR #151 at fd22f71. This cumulative PR targets main and remains unmerged with the feature stack.

Write one tagged HDF5 response file for every stage-2 reconstructed propagator. Each file retains the co-passing and signed counter-passing current-row contractions before spatial quadrature, the separate directional weights, geometry, passing-grid size, and the weighted local qflux partition.

The globally reconstructed incoming boundary rows remain authoritative. Summing all 199 tagged partitions reproduces the final joined current-row qflux, including the response represented by join_ripples and join_ends.

Physics and numerical invariants:

  • The kinetic operator, source, Runge-Kutta stencil, eta grid, spatial quadrature, join algebra, convergence criteria, ABI, and memory layout are unchanged.
  • The current row retains the existing co-minus-counter sign and solver force order A1,A3,A2.
  • Boozer coordinates, CGS normalization, Fourier phase, boundary conditions, and final gamma_out/D31/D32 are unchanged.
  • The output is a reconstructed normalized response coefficient. It is not yet physical pointwise current: adjacent one-sided interface values differ by 2.5% and 5.7% of the global A1/A2 amplitudes, and physical forces have not been applied.

Verification

Test fails on main

$ find cases/nminus_phi -name 'local_parallel_response_*.h5' | wc -l
0
$ python3 analyze_local_response.py cases/nminus_phi results
FileNotFoundError: no local response files

The previous single-call qflux_current_contributions.dat is replaced on every ripple call and cannot decompose final joined qflux.

Test passes after fix

$ fo
Static: OK (343 modules)
Build: OK
Tests: OK
Lint: OK
All stages passed

$ python3 run_reconstruction.py /path/to/NEO-2-local-parallel-response
propagator_files: 398
reconstructed_boundaries: 199
local_response_files: 199
runtime_markers: 0

qflux(2,1) relative closure error: 4.806611157179914e-15
qflux(2,3) relative closure error: 8.636816603296595e-15
point rows: 99415

The focused qflux_profile test passes. All added lines are at most 88 columns; the two newly formatted files pass fo fmt --check.

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.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant