Repo - API! - Release GPEC v2.0.0

.JuliaFormatter.toml

@@ -3,10 +3,10 @@ always_for_in = true
 for_in_replacement = "in"
 annotate_untyped_fields_with_any = false
 whitespace_in_kwargs = false
-always_use_return = true
+always_use_return = false
 format_docstrings = true
 remove_extra_newlines = false
 trailing_comma = false
 join_lines_based_on_source = true
 separate_kwargs_with_semicolon = true
-whitespace_ops_in_indices = false
+whitespace_ops_in_indices = false

.claude/agent-memory/fortran-physics-reviewer/MEMORY.md

@@ -0,0 +1,8 @@
+# Agent Memory Index
+
+- [Fortran-Julia Correspondence Map](fortran_correspondence_map.md) — Maps Fortran GPEC source files to Julia JPEC modules
+- [Known Physics Issues in PerturbedEquilibrium](known_issues_perteq.md) — xi_modes placeholders, the xi_n Jacobian shortcut
+- [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify
+- [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify
+- [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS
+- [reg_spot vs singfac_min audit](reg_spot_regularization.md) — reg_spot field-recon smoothing CORRECT; distinct from singfac_min ODE gate

.claude/agent-memory/fortran-physics-reviewer/feedback_annotations.md

@@ -0,0 +1,13 @@
+---
+name: annotation_style
+description: How to cite Fortran code in annotations — use papers and submodule names, not line numbers
+type: feedback
+---
+
+Do NOT cite Fortran source line numbers in code annotations (e.g., "gpeq.f:100"). The Fortran codebase is actively developed and line numbers change. Instead, cite:
+- Paper equations (e.g., "Park et al. 2007, Eq. 8")
+- Fortran subroutine/function names (e.g., "matches Fortran `gpeq_contra`")
+
+**Why:** Fortran GPEC is still under active development; line-number references rot quickly.
+
+**How to apply:** When the fortran-physics-reviewer suggests line-number annotations, convert them to subroutine-name or paper-equation references before writing to code.

.claude/agent-memory/fortran-physics-reviewer/fortran_correspondence_map.md

@@ -0,0 +1,43 @@
+---
+name: Fortran-Julia File Correspondence Map
+description: Mapping between Fortran GPEC source files and Julia JPEC modules for physics review
+type: reference
+---
+
+Fortran paths are relative to the root of the Fortran GPEC clone; Julia paths are
+relative to this repository.
+
+## Coil/ForcingTerms
+- `coil/coil.F` (coil_read) -> `src/ForcingTerms/CoilGeometry.jl` (apply_transforms)
+- `coil/field.F` (field_bs_psi) -> `src/ForcingTerms/BiotSavart.jl` + `src/ForcingTerms/CoilFourier.jl`
+
+## PerturbedEquilibrium
+- `gpec/gpeq.f` (gpeq_sol, gpeq_contra, gpeq_surface, gpeq_normal) -> `src/PerturbedEquilibrium/FieldReconstruction.jl` + `src/PerturbedEquilibrium/ResponseMatrices.jl`
+- `gpec/gpresp.f` (gpresp_pinduct, gpresp_sinduct, gpresp_permeab) -> `src/PerturbedEquilibrium/ResponseMatrices.jl`
+- `gpec/gpout.f` (gpout_singcoup, gpout_xbnormal) -> `src/PerturbedEquilibrium/SingularCoupling.jl` + `src/PerturbedEquilibrium/FieldReconstruction.jl`
+- `gpec/gpvacuum.f` (gpvacuum_flxsurf) -> `src/Vacuum/Vacuum.jl` (compute_Iv branch of _compute_vacuum_response_2d!) + `src/PerturbedEquilibrium/ResponseMatrices.jl` (calc_surface_inductance)
+
+## KineticForces (NTV, Fortran `pentrc/`) — see kinetic_ntv_map.md for the audit checklist
+- `pentrc/torque.F90` -> `src/KineticForces/Torque.jl` (tpsi! single-surface torque)
+- `pentrc/pitch.f90` -> `src/KineticForces/PitchIntegration.jl` (pitch-angle lambda integration, bounce-averaged operator assembly)
+- `pentrc/energy.f90` -> `src/KineticForces/EnergyIntegration.jl` (energy-space quadrature integrand)
+- `pentrc/pentrc.F90` (orchestration) -> `src/KineticForces/{KineticForces.jl,Compute.jl}`
+- `pentrc/inputs.f90` + `params.f90` -> `src/KineticForces/KineticForcesStructs.jl`
+- `pentrc/dcon_interface.f` -> `src/KineticForces/CalculatedKineticMatrices.jl` (bridge to ForceFreeStates kinetic stability) + `src/ForceFreeStates/{Kinetic.jl,FixedKineticMatrices.jl}`
+- bounce averaging of the 6 perturbed-action matrices -> `src/KineticForces/BounceAveraging.jl`
+
+## InnerLayer (resistive matched-asymptotics, Fortran `rmatch/`) — see resistive_layer_map.md for the audit checklist
+- `rmatch/deltac.f` (Galerkin solver) -> `src/InnerLayer/GGJ/Galerkin.jl`
+- `rmatch/deltar.f` (shooting solver) -> `src/InnerLayer/GGJ/Shooting.jl`
+- `rmatch/inps.f` + `inpso.f` (Wasow asymptotic basis: T,J,P,B,Q,C,D,Y,Z,U matrices) -> `src/InnerLayer/GGJ/InnerAsymptotics.jl`
+- `rmatch/{inner.f,match.f,msing.f,gamma.f}` -> `src/InnerLayer/GGJ/` (parameters, matching data, special functions)
+- SLAYER drift-MHD two-fluid solver: `src/InnerLayer/SLAYER/` is a placeholder pending implementation
+
+## Key Fortran Conventions
+- `sq%f(1)` = F (toroidal field function = R*B_tor), `sq%f(4)` = q (safety factor)
+- `sq%f1(4)` = dq/dpsi
+- `chi1 = twopi*psio` (poloidal flux normalization)
+- `singfac = mfac - nn*q` (resonance factor)
+- `ifac = (0,1)` (imaginary unit)
+- `wegt=0`: no weighting; `wegt=1`: J*|grad psi| weighting (default for bmn)
+- Vacuum Green's functions use reversed theta in Fortran (rtheta = mthsurf - itheta)

.claude/agent-memory/fortran-physics-reviewer/galerkin_assembly_map.md

@@ -0,0 +1,83 @@
+---
+name: galerkin-assembly-map
+description: RDCON singular Galerkin Δ′ solver — gal.f ↔ GalerkinAssembly.jl element-level weak-form correspondence and audit notes
+metadata:
+  type: project
+---
+
+# RDCON Galerkin Δ′ assembly: gal.f ↔ GalerkinAssembly.jl
+
+File: `src/ForceFreeStates/Galerkin/GalerkinAssembly.jl` ↔ `GPEC_fortran/rdcon/gal.f`.
+On branch `feature/rdcon-galerkin`. Reviewed 2026-06-16.
+
+## Correspondence (verified, not just docstring-trusted)
+- gal_hermite (222-262) ↔ gal_hermite: pb/qb EXACT.
+- gal_get_fkg (1013-1071) ↔ gal_get_fkg: F=Fbar·sf[i]·sf[j], K=K·sf[i], G=G. singfac=mlow-nn*q+(0..mpert-1) = DIRECT (m-nq), not reciprocal. Correct.
+- gal_gauss_quad (1079-1138) ↔ gal_gauss_quad!: F·qq + K·qp + conj(K[jpert,ipert])·pq + G·pp. swap pb(2)/pb(3) (0-based) = Julia pb[3]/pb[4]. Correct.
+- gal_extension (1146-1322) ↔ gal_extension!: emat/ediag/rhs/erhs all match incl surface terms 1277-1300 and w1/w2.
+- gal_resonant! ↔ gal_lsode_int(780-900)+gal_lsode_der(908-954)+gal_make_arrays signs(1366-1369). QuadGK replaces LSODE.
+- gal_assemble_mat (694-772) ↔ gal_assemble_mat!: LU offset=kl+ku+1, chol offset=1, conj(e) Hermitian placement. Correct.
+- gal_assemble_rhs (642-686) ↔ gal_assemble_rhs!: isol++ at ext2-left / res-right. Correct.
+- gal_set_boundary (1554-1620) ↔ gal_set_boundary!: axis idx, edge 3-way (rpec/vac/fixed). DOF index 3(0-based)=Julia 4.
+
+## Resonant sign chain (the trap — verified correct)
+Fortran: LSODE integrates x0→x_lsode giving raw u; then line 882 `IF right: u=-u`; then make_arrays:
+erhs=-u_res(1), ediag=+u_res(2), rhs=-u_hermite1, emat=+u_hermite2.
+Net with s=+1(left)/-1(right): erhs=-s·res1, ediag=s·res2, rhs=-s·hbig, emat=s·hsmall.
+Julia sgn_u = s (=-1 right/+1 left): erhs=-sgn_u·raw1, ediag=sgn_u·raw2, rhs=-sgn_u·hbig, emat=sgn_u·hsmall. MATCHES.
+QuadGK limits x0→x_lsode preserve LSODE orientation so raw == pre-negation u. Correct.
+du_res conj-weight uses the SMALL solution (isol(2)) in both. Correct.
+
+## Only open item (MINOR, out-of-file)
+sing_get_ua_gal/dua_gal two-sided z<0 → sig=-1 series at -z with derivative ×(-1) is a REIMPLEMENTATION
+of sing.f's internal singular-surface handling, not a line port. The ×(-1) chain is internally consistent
+but worth a dedicated check against sing.f sing_get_ua/dua if Δ′ ever disagrees with Fortran. Lives outside
+GalerkinAssembly.jl.
+
+## RPEC matching (GalerkinMatch.jl <-> rmatch/match.f match_rpec) — reviewed 2026-06-20
+- match block (mat idx1..idx4, -delta1/delta2 signs) + rmat=-TRANSPOSE(delta coil block) + cout/cin split: EXACT port. PASS.
+- matched outer soln Σcout(isol)·sols + sols(:,:,csol): EXACT (match_output_solution L32-41).
+- resist.f E/F/H/M/G/K/taua/taur/v1: EXACT (Resist.jl). eta/rho caller-supplied = deliberate documented deviation.
+- gal quadrature: jacobi_alloc "gll" = Gauss-Lobatto = Julia gausslobatto. nq1 lower=|nn·q1|, upper=|q1|: match.
+
+## SUSPECT (should-fix) — RPEC forced eigenvalue γ
+GalerkinMatch.jl:66 uses γ = 2π·im·n·gal_rotation. Fortran match_rpec passes `initguess` (namelist,
+default 0.0) DIRECTLY to deltac_run — NOT i·n·rotation. The `guess+ifac*ntor*rotation` Doppler form
+lives only in match_delta (Newton eigenvalue search), NOT the RPEC path. Two issues:
+  (1) extra 2π: Fortran eig=ifac*ntor*rotation has NO 2π (rotation is ω[rad/s], eig is 1/time, q0=x0/taua).
+  (2) RPEC path in Fortran uses initguess (the user forced eigenvalue), not a per-surface i·n·rotation shift.
+Confirm intended convention vs Fortran before trusting resistive Δ(Q); ideal_flag path (γ unused) is unaffected.
+RESOLUTION 2026-06-20: author confirms the Hz convention is DELIBERATE (gal_rotation is f[Hz], so γ=2πi·n·f);
+the overstated comments were corrected. Δ(Q) now guarded by test/runtests_innerlayer.jl (Julia self-pin on
+glasser_wang_2020_eq55). Remaining follow-up: upgrade that guard to a Fortran deltac_run / paper
+cross-check to fully close the convention question.
+
+## Verdict: PASS WITH REQUIRED ANNOTATIONS for outer solve; RPEC γ convention is SHOULD-FIX (verify vs match.f).
+
+
+# Second pass
+
+
+## PR #266 "Galerkin Integrator" full review (2026-06-25) — PASS
+Exhaustive line-by-line vs gal.f / match.f / deltac.f / resist.f / sing.f. All physics-bearing code faithful.
+- RESOLVED (was SUSPECT): forced eig γ=2πi·n·rotation IS correct. match.f match_rpec line:
+  `rpec_eigenvalues(ising)=2*pi*rotation(ising)*REAL(ntor,r8)*ifac` — the 2π IS in Fortran. Julia matches exactly.
+  My earlier "no 2π in Fortran" note was wrong (confused match_delta Newton path with match_rpec). Closed.
+- Verified EXACT: gal_get_fkg (F=QF̄Q,K=QK̄,G=Ḡ,singfac=m-nq), gauss_quad weak form (F qq+K qp+K† pq+G pp),
+  gal_extension, gal_resonant (mpert-fastest hermite layout; Fortran's apparent (0:np,mpert,2) reshape is a
+  no-op relabel — der writes mpert-fastest, caller reads mpert-fastest), assemble_mat/rhs, set_boundary,
+  gal_get_solution (restore_uh/us/ul), PEST-3 A/B/Γ/Δ ± combos, Δ′ extraction, match_rpec matrix signs,
+  resist_eval E/F/G/H/K/M/taua/taur, mercier_di=E+F+H-1/4, rescale=sfac^(2p1/3)·v1^(2p1) + delta(1)↔(2) swap
+  (swap applied in BOTH solve_inner and solve_inner_profile), sing_matvec, sing_get_dua, ξ_s=-A⁻¹(Bξ'+Cξ).
+- Documented deviations (all fine): η/ρ caller-supplied (no Spitzer); gal_gamma default 5/3 = resist.f;
+  QuadGK replaces LSODE; sing_get_dua omits Fortran `sig` (moved to sing_get_dua_gal left-side ×(-1)).
+- InnerAsymptotics.jl + Shooting.jl diffs = PURE JuliaFormatter whitespace, zero numeric change.
+- Reference.jl G=8.950e1,H=1.292e-2,K=2.332e2 = corrected GW2020 Eq.55 (prior transcription bug fixed).
+
+## OPEN ITEM (latent, n≥2 only): gal_make_grid resonant-cell width asymmetry
+gal.f gal_make_grid genuinely uses nq1=ABS(nn*sing(ising)%q1) on the LOWER bound but nq1=ABS(sing(ising+1)%q1)
+(NO nn) on the UPPER bound. Julia GalerkinGrid.jl faithfully reproduces. Natural layer scale is 1/(n q1) on BOTH
+sides (singfac'≈-n q1), so the upper/left side of each interior surface gets cells ~n× too wide for n≥2.
+No effect at n=1 (RDCON resistive target). galerkin_solve allows single n≥2 (blocks only multi-n npert>1), so it
+IS reachable. Recommend: keep faithful but warn, or fix both sides to nn*q1 (deviation, more correct) gated by
+validation. Suspected upstream Fortran oversight, not a Julia port error.

.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md

@@ -0,0 +1,90 @@
+---
+name: KineticForces (NTV) Audit Checklist
+description: Curated map and physics-audit checklist for the KineticForces module (NTV, formerly PENTRC); what to read and what to verify
+type: reference
+---
+
+The KineticForces module computes neoclassical toroidal viscosity (NTV) torque and the
+kinetic energy/matrices from trapped-particle nonambipolar transport. Fortran counterpart is
+the Fortran `pentrc/` sources (file mapping in fortran_correspondence_map.md).
+
+## Governing theory
+- Logan & Park (2013) PoP 20, 122507: diamagnetic frequency (Eq. 7), energy integrand (Eq. 8), torque normalization Im(T) = 2n·δW_k (Eq. 19).
+- Logan (2015) PhD Thesis, Ch. 7: the six kinetic matrices A_k, B_k, C_k, D_k, E_k, H_k (Eqs 7.30–7.35) as energy-space integrals of perturbed action operators W_X, W_Y, W_Z; resonance splitting/suppression F_h = (Q−P†)F̄(Q−P)+… (Eq. 7.46). App. C: DCON matrix form of the perturbed action. App. D: numerical treatment of integrable singularities in bounce averages.
+- Park et al. (2009) PRL 102, 065002: trapped-particle nonambipolar transport foundation.
+
+## Key Julia files and roles
+- `Torque.jl` — `tpsi!()`: single-surface torque; poloidal grid sampling; diamagnetic frequencies (Eq. 7).
+- `EnergyIntegration.jl` — energy-space quadrature integrand (Eq. 8); Maxwellian vs JKP distribution options.
+- `PitchIntegration.jl` — pitch-angle (λ) integration; bounce-averaged operator assembly.
+- `BounceAveraging.jl` — bounce-averaging of the perturbed-action matrices (Eqs 7.30–7.35); packs 3 Hermitian + 3 full blocks.
+- `Compute.jl` — `compute_torque_all_methods!()` orchestration over all ψ surfaces.
+- `CalculatedKineticMatrices.jl` — bridge to ForceFreeStates kinetic stability (`Kinetic.jl`, `FixedKineticMatrices.jl`).
+- `KineticForcesStructs.jl` — `KineticForcesControl` (TOML params) and `KineticForcesInternal` (state, hints, mode indexing).
+- `Output.jl` — HDF5 writer.
+
+## What to verify in a review
+- **Matrix structure**: A symmetric; the Hermitian blocks actually Hermitian; the 3-Hermitian + 3-full packing in BounceAveraging matches Logan 2015 Eqs 7.30–7.35.
+- **Resonance handling**: Sokhotski–Plemelj pole shift / singular-denominator gating near rational surfaces (App. D). Check the singular-eps threshold is applied where the denominator vanishes, not globally.
+- **Quadrature tolerances**: outer ψ (atol_psi/rtol_psi) and inner pitch/energy (atol_xlmda/rtol_xlmda) are passed through, not hard-coded to lazy defaults.
+- **Distribution/collision options**: energy integrand (Eq. 8) honors the configured collision operator (harmonic/Krook/zero) and distribution (Maxwellian/JKP/CGL) — not silently fixed to one.
+- **Normalization**: torque Im(T)=2n·δW_k (Eq. 19) and diamagnetic-frequency sign/factor conventions match Fortran `torque.F90`.
+- **Mode indexing**: m, n ranges and block packing over n stay consistent with ForceFreeStates.
+- **Method variants**: FGAR/TGAR/PGAR/RLAR/CLAR/FCGL/TMM/WMM each present, not stubbed to a single fallback.
+
+## Multi-ion (D-T) NTV — composition/collisionality (2026-07)
+Both Fortran PENTRC (`inputs.f90:236-243`, `read_kin`) and Julia (`KineticProfiles.jl:261-278`)
+support ONE main ion (zi,mi) + ONE impurity (zimp,mimp) per run. Correct multi-main-ion NTV is
+an EXTENSION beyond Fortran, but consistent with Logan-Park 2013's pitch-angle (Lorentz) model.
+- Zeff = Σ_s Z_s² n_s / n_e; quasineutrality n_e = Σ_s Z_s n_s (all ions incl. impurity).
+- Bug in current D-T split (run D, then T, each with ni=Ni/2): the `ni` column sets BOTH species
+  density AND Zeff via z = zimp-(n_i/n_e)zi(zimp-zi). With ni=Ni/2 it treats the missing half as
+  high-Z impurity → Zeff≈3.7 instead of true ≈1–1.4. Corrupts zpitch → nue,nui.
+- Correct collisionality: ν_a ∝ Z_a² lnΛ · Σ_b n_b Z_b² /(√m_a T_a^{3/2}) = ∝ Z_a² lnΛ·(n_e·Zeff_true).
+  The zpitch(Zeff) polynomial is a main-ion+impurity closure (momentum-restoring correction); the
+  MINIMAL fix is to feed the TRUE (full-composition) Zeff into the existing zpitch/ν formulas.
+- Additivity: τ = Σ_s τ_s. Lorentz operator is additive over field species; species couple only
+  through shared δB, shared ω_E, shared Zeff/ν. Additive at this theory's order.
+- What changes single→multi: Zeff, zpitch, nue, nui (→nueff). UNCHANGED for equal-shape D-T:
+  wdian/wdiat (log-derivative, density factor cancels), wtran/wbhat/wdhat/wgyro (already per-species),
+  and the resonant-density prefactor (Ni/2 per species is correct).
+- ASIDE (separate fidelity bug, same code block): Fortran `inputs.f90:238` uses natural log for lnΛ;
+  Julia `KineticProfiles.jl:270` uses log10 — diverges away from the n=1e20,T=1keV reference point.
+
+## Multi-ion NTV — full-composition species set (2026-07 audit, PASS-with-caveats)
+Reviewed `resolve_ntv_species` (KineticProfiles.jl ~214-266) + `compute_calculated_kinetic_matrices`
+(CalculatedKineticMatrices.jl ~96-151). Verdict: physics is sound.
+- ν_s field-density RECONCILIATION (supersedes the earlier "should be n_e·Zeff" note above):
+  code uses `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ/(√m_s·T_i^1.5)`, i.e. field density = zpitch·n_main
+  (n_main = Σ MAIN-ion densities, no impurity, unweighted). This is CORRECT and MORE faithful to
+  single-ion PENTRC than n_e·Zeff: PENTRC's design is zpitch·n_i, NOT Σ_b n_b Z_b². The two are NOT
+  numerically equal (Zeff=1.5, C6, n_i/n_e=0.9 → zpitch·n_i = 1.225·n_e vs n_e·Zeff = 1.5·n_e, ~18%
+  apart); that gap IS the intended momentum-restoring design difference that zpitch(Zeff) carries.
+  n_main vs n_e·Zeff is a deliberate PENTRC-fidelity choice, exact in the z=1 main-ion domain.
+- z_s² test-particle factor: CORRECT and correctly placed (deflection freq ∝ test charge²). Single-ion
+  had no z² only because zi=1. Reduces EXACTLY to single-ion nui for one z=1,fraction=1 ion (verified).
+- Impurity as its own test species: field density zpitch·n_main is NOT undercounting — the impurity's
+  z_imp² contribution is already folded into Zeff inside zpitch. CAVEAT: zpitch is strictly a main-ion
+  momentum-restoring closure; reusing it for the impurity/electron ν is an approximation beyond the
+  single-ion theory. Acceptable (impurity δW ∝ n_imp is small); worth a one-line annotation.
+- Electron descriptor passes `ns[1]` (first ion's density) as its `ni_spline` — HARMLESS: the kernel
+  `_setup_surface_state` (Torque.jl:659-664) reads `ne_spline` (full shared n_e) when electron=true,
+  never ni_spline. Electron n_s = full n_e (correct). Cosmetic smell only.
+- Self-consistent δW summation (CalculatedKineticMatrices.jl:112-148): kw_flat/kt_flat are PURELY the
+  kinetic matrices (Logan 7.30-7.35), every term ∝ species phase-space density n_s·f0_s — NO species-
+  independent baseline. Fluid F,K,G added ONCE downstream in _compute_fkg_matrices!, outside the
+  species loop. So `+=` over species is clean additivity; NOTHING in kw is wrongly ×species-count.
+  → TC-24 n=3 δW +0.066(D) → −0.10(D+T+C+e) sign flip is PLAUSIBLE physics near marginal stability,
+  NOT a double-count. Dominant driver: the newly-ON electron channel (full n_e, opposite precession).
+  DECISIVE cheap diagnostic to confirm: run D(½)+T(½) with electron=OFF, impurity absent — should
+  return ≈ +0.066 (single-ion D). If D+T alone ≈ +0.066, the whole shift is electron+impurity = physical.
+
+## Single-ion nui vs multi-species z_s² (#339, 2026-07 decision)
+Fortran PENTRC `inputs.f90:240-241` single-ion nui = (zpitch/3.5e17)·n_i·lnΛ/(√mi·T_i^1.5) has
+NO explicit zi² (implicitly assumes zi=1, main ion hydrogenic). Julia `load_kinetic_profiles`
+(~L416) is a faithful exact port. Multi-species `_nu` (~L253) adds explicit test-particle z_s²
+(the physically correct pitch-angle/Lorentz form; zpitch is the field-side momentum-restoring
+factor of Zeff, independent of test charge — no double count). RECOMMENDATION: option (a) ADD
+zi² to single-ion nui. It is numerically identical for the default zi=1 (regression byte-identical),
+makes single-ion the true 1-species special case of `_nu`, and only "deviates" from Fortran in the
+exotic zi≠1 case where Fortran is physically wrong anyway. Annotate as a documented improvement.