KineticForces Module (formerly PENTRC)

CLAUDE.md

@@ -95,9 +95,9 @@ GPEC will eventually implement resistive MHD stability analysis based on:
   - Published: Physics of Plasmas **27**, 122509 (2020)
   - Describes: Asymptotic matching for resistive plasma response
 
-### PENTRC Module (Future Work)
+### KineticForces Module (NTV)
 
-GPEC will eventually port the PENTRC (Perturbed Equilibrium Neoclassical Toroidal viscosity in Realistic geometry Code) functionality from the Fortran GPEC suite. This is described in:
+The KineticForces module (formerly PENTRC) implements neoclassical toroidal viscosity calculations. Based on:
 
 - **Logan & Park (2013)**: "Neoclassical toroidal viscosity in perturbed equilibria with general tokamak geometry"
   - Location: `docs/resources/2013-Logan-Neoclassical_toroidal_viscosity_in_perturbed_equilibria_with_general_tokamak_geometry.pdf`
@@ -107,7 +107,7 @@ GPEC will eventually port the PENTRC (Perturbed Equilibrium Neoclassical Toroida
 - **Logan (2015)**: "Electromagnetic Torque in Tokamaks with Toroidal Asymmetries"
   - Location: `docs/resources/2015-Logan-Electromagnetic_Torque_in_Tokamaks_with_Toroidal_Asymmetries-compressed.pdf`
   - Published: PhD Thesis, Princeton University (2015)
-  - Describes: Complete PENTRC theory and implementation. **Chapter 7** details the hybrid drift-kinetic MHD eigenfunction calculation: 6 kinetic matrices Ak,Bk,Ck,Dk,Ek,Hk (Eqs 7.30-7.35) as energy-space integrals of perturbed action operators WX,WY,WZ; hybrid Euler-Lagrange equations; resonance splitting/suppression where Fh=(Q-P†)F̄(Q-P)+... shifts singularities away from rational surfaces (Eq 7.46); convergence to ideal limit. **Appendix C** derives the DCON matrix form of the perturbed action (Eqs C.1-C.11) used to compute the kinetic coefficient matrices. **Appendix D** details numerical treatment of integrable singularities in bounce averages.
+  - Describes: Complete NTV theory and implementation. **Chapter 7** details the hybrid drift-kinetic MHD eigenfunction calculation: 6 kinetic matrices Ak,Bk,Ck,Dk,Ek,Hk (Eqs 7.30-7.35) as energy-space integrals of perturbed action operators WX,WY,WZ; hybrid Euler-Lagrange equations; resonance splitting/suppression where Fh=(Q-P†)F̄(Q-P)+... shifts singularities away from rational surfaces (Eq 7.46); convergence to ideal limit. **Appendix C** derives the DCON matrix form of the perturbed action (Eqs C.1-C.11) used to compute the kinetic coefficient matrices. **Appendix D** details numerical treatment of integrable singularities in bounce averages.
 
 ### Additional References
 
@@ -584,6 +584,7 @@ This format is used for compiling release notes, so tags should be human-readabl
 - **Indexing**: The codebase uses 0-based indexing in many places to match Fortran conventions, then converts to 1-based Julia indexing
 - **No step numbering in code comments** - Avoid annotations like "Step 1: do this" followed by "Step 2: do that". These get out of sync as code changes. Just describe the action without numbering.
 - **Documentation coverage** - When adding a new module or submodule with public docstrings, add a corresponding `@autodocs` block in `docs/src/`. Documenter CI will fail with a `missing_docs` error if any exported docstring is not covered. The analysis submodule docs live in `docs/src/analysis.md`.
+- **Docstrings are rendered as Markdown** - Documenter parses docstrings as CommonMark, so `[text](...)` patterns become hyperlinks and will fail CI with `invalid local link/image` if the target doesn't exist. Common pitfall: unit annotations like `[degrees] (or [m] if ...)` parse as `[degrees](or [m] if ...)` — a broken link. Use plain words (`in degrees`) or backticks (`` `[m]` ``) for unit labels inside docstrings. Same rule for bracketed array-shape hints (`[ncoil]`) followed by parentheses. When in doubt, preview with `julia --project=docs docs/make.jl` before pushing.
 - **Keep code comments concise** - A comment should be one line where possible. Do not write multi-line block comments explaining the current session's investigation, what was tried, what was wrong before, or why a specific file/path behaves differently. State what the code does and why at a general level. Example of too much detail: a 6-line block explaining that efit_by_inversion uses psilow>0 while CHEASE starts at 0, that the old code was removed, and that spline spikes result. Preferred: `# Replicate Fortran inverse.f: overwrite deta at axis (r²=0) by extrapolating from innermost surfaces.`
 
 ### Output Files
@@ -607,6 +608,19 @@ This format is used for compiling release notes, so tags should be human-readabl
   plot!(p, m_ext, a_ext; seriestype=:steppre, lw=2, label="...")
   ```
 
+### Subagent Consultations
+
+When delegating to specialized agents (julia-performance-optimizer, fast-interpolations-optimizer, fortran-physics-reviewer, clean-code-reviewer, etc.), every prompt **must include an explicit budget** because runaway agents silently consume the user's daily token quota. A single consultation that explores instead of editing has been measured at 167 tool calls / 55 minutes / no return — that is a session-killer. Defaults:
+
+- **Hard cap: ≤ 30 tool uses and ≤ 10 minutes wall time per consultation.** State both numbers in the prompt verbatim ("Budget: ≤30 tool uses, ≤10 min").
+- **Single concrete deliverable.** One file, one function, or one named hotspot list. Not "audit the module."
+- **No exploration phase.** The prompt must hand the agent the file paths and line numbers; the agent's job is to edit, not to map the codebase.
+- **Require an interim status if the work might exceed budget.** Tell the agent: "If you cannot finish within budget, stop and report what was changed and what remains."
+- **Prefer two short focused agents over one open-ended one.** If an investigation needs both performance and interpolation review, run them sequentially with separate ≤30-tool budgets — don't chain them in one long prompt.
+- **Never re-launch a runaway agent.** If an agent hits the API rate limit before returning, do not retry; report the partial state to the user and switch to hand-implementation.
+
+These rules apply to **every** Agent tool invocation, not just performance work.
+
 ### Code Formatting
 
 Pre-commit hooks enforce formatting via JuliaFormatter (v1.0.62) and general file hygiene. **All code you write or modify must already conform to these standards before committing**, so the hooks have nothing to fix. Failing to do this creates noisy diffs in PRs where formatting changes leak into unrelated files.

benchmarks/_plot_cond_fbar.jl

@@ -0,0 +1,28 @@
+"""
+    _plot_cond_fbar.jl
+
+Thin shim used by the kinetic-DCON benchmark scripts. Delegates to
+`GeneralizedPerturbedEquilibrium.Analysis.ForceFreeStates.plot_cond_fbar`
+so the benchmarks share the same visualisation that every user gets from
+the Analysis module.
+"""
+
+using GeneralizedPerturbedEquilibrium
+const _AnalysisFFS = GeneralizedPerturbedEquilibrium.Analysis.ForceFreeStates
+
+"""
+    plot_cond_fbar_scan(h5path, png_path; title="") → png_path or nothing
+
+Save a cond(F̄) vs ψ plot to `png_path` by calling
+`Analysis.ForceFreeStates.plot_cond_fbar(h5path; save_path=png_path)`. The
+`title` argument is ignored (the Analysis function composes its own title
+from the stored kmsing count); the shim keeps the old signature so
+existing benchmark scripts do not need to change.
+"""
+function plot_cond_fbar_scan(h5path::String, png_path::String; title::String="")
+    _ = title  # kept for backward compatibility; Analysis function composes its own title
+    p = _AnalysisFFS.plot_cond_fbar(h5path; save_path=png_path)
+    p === nothing && return nothing
+    println(stderr, "  cond(F̄) plot: ", abspath(png_path))
+    return abspath(png_path)
+end

benchmarks/benchmark_a10_kinetic_dcon.jl

@@ -0,0 +1,298 @@
+#!/usr/bin/env julia
+"""
+    benchmark_a10_kinetic_dcon.jl
+
+Fast-iteration kinetic-DCON benchmark against the Fortran GPEC
+`a10_kinetic_example` (mpsi=16 runs in ≈48 s). Mirrors the structure of
+`benchmark_diiid_kinetic_dcon.jl` but at the smaller a10 problem size so
+each iteration of the KF→FFS self-consistent bug-hunt stays under ~10 min
+of wall time.
+
+Compares (least-stable mode, index 1):
+  - Re(et[1]) relative error
+  - Im(et[1]) relative error
+  - |⟨ξ_J | ξ_F⟩| total-energy eigenvector overlap
+
+Inputs are read directly from the user's local Fortran GPEC checkout; no
+a10 inputs are duplicated into this repo.
+
+Acceptance (from validated-exploring-clarke plan):
+  - |Re err| < 10 %
+  - |Im err| < 30 %
+  - |⟨ξ_J | ξ_F⟩| > 0.90
+
+Usage:
+    julia --project=. benchmarks/benchmark_a10_kinetic_dcon.jl [fortran_kinetic_dir]
+
+    fortran_kinetic_dir defaults to \$GPEC_FORTRAN_A10_DCON or
+    ~/Code/gpec/docs/examples/a10_kinetic_example.
+"""
+
+using Printf
+using HDF5
+using NCDatasets
+using LinearAlgebra
+
+using GeneralizedPerturbedEquilibrium
+const GPE = GeneralizedPerturbedEquilibrium
+
+include(joinpath(@__DIR__, "_plot_cond_fbar.jl"))
+
+const DEFAULT_FORTRAN_KIN_DIR   = expanduser("~/Code/gpec/docs/examples/a10_kinetic_example")
+const DEFAULT_FORTRAN_IDEAL_DIR = expanduser("~/Code/gpec/docs/examples/a10_ideal_example")
+
+"""
+    discover_inputs(fortran_kin_dir) → (; eq_file, kin_file, dcon_nc)
+
+Locate the EFIT g-file and kinetic profile for the a10 example (both live
+in `a10_ideal_example`), plus the kinetic DCON NetCDF reference from
+`a10_kinetic_example`. Errors if any required file is missing.
+"""
+function discover_inputs(fortran_kin_dir::String)
+    isdir(fortran_kin_dir) || error("Fortran a10 kinetic example directory not found: $fortran_kin_dir")
+    # a10 shares its g-file and .kin with the ideal example directory.
+    ideal_dir = normpath(joinpath(fortran_kin_dir, "..", "a10_ideal_example"))
+    isdir(ideal_dir) || error("a10 ideal example directory not found: $ideal_dir")
+    eq_file  = joinpath(ideal_dir, "fix_a100_k10_q2_bn010_prof1")
+    kin_file = joinpath(ideal_dir, "a10_prof1.txt")
+    dcon_nc  = joinpath(fortran_kin_dir, "dcon_output_n1.nc")
+    for (label, path) in [("EFIT g-file", eq_file),
+                          ("kinetic profile", kin_file),
+                          ("dcon NetCDF", dcon_nc)]
+        isfile(path) || error("$label not found: $path")
+    end
+    return (; eq_file, kin_file, dcon_nc)
+end
+
+"""
+    read_dcon_reference(dcon_nc) → NamedTuple
+
+Read least-stable (mode=1) total-energy eigenvalue and eigenvector from
+the Fortran a10 `dcon_output_n1.nc`. NCDatasets reverses NetCDF dim order
+so Fortran `W_t_eigenvalue(i, mode)` → Julia `Wt[mode, i]` and
+`W_t_eigenvector(i, mode, m)` → Julia `Wte[m, mode, i]`. `i ∈ {1=Re, 2=Im}`.
+Global attributes `mlow`, `mhigh`, `mpert` give the Fortran m-grid so the
+Julia-side eigenvector can be aligned by m index.
+"""
+function read_dcon_reference(dcon_nc::String)
+    NCDatasets.Dataset(dcon_nc, "r") do ds
+        Wt  = Array(ds["W_t_eigenvalue"][:, :])        # Julia (mode, i)
+        Wte = Array(ds["W_t_eigenvector"][:, :, :])    # Julia (m, mode, i)
+        Wp  = Array(ds["W_p_eigenvalue"][:, :])
+        Wv  = Array(ds["W_v_eigenvalue"][:, :])
+        mlow  = Int(ds.attrib["mlow"])
+        mhigh = Int(ds.attrib["mhigh"])
+        mpert = Int(ds.attrib["mpert"])
+        # Least-stable eigenvector: all m, mode=1, re+im → length mpert
+        xi_F = complex.(Wte[:, 1, 1], Wte[:, 1, 2])
+        return (
+            W_t_re = Wt[1, 1], W_t_im = Wt[1, 2],
+            W_p_re = Wp[1, 1], W_p_im = Wp[1, 2],
+            W_v_re = Wv[1, 1], W_v_im = Wv[1, 2],
+            xi_F   = xi_F,
+            mlow   = mlow,
+            mhigh  = mhigh,
+            mpert  = mpert,
+        )
+    end
+end
+
+"""
+    build_benchmark_tomldir(eq_file, kin_file) → tmpdir
+
+Construct a temporary directory containing a Julia `gpec.toml` with a10
+kinetic-DCON settings mirroring `a10_kinetic_example/{dcon.in, equil.in,
+pentrc.in}`, plus a symlink (or copy) to the EFIT g-file. The
+`[KineticForces]` section points at the Fortran `.kin` file by absolute
+path. The term-by-term mapping is also committed to
+`examples/a10_kinetic_example/gpec.toml`.
+"""
+function build_benchmark_tomldir(eq_file::String, kin_file::String)
+    tmpdir = mktempdir(; prefix="a10_kinetic_dcon_bench_")
+    eq_name = basename(eq_file)
+    try
+        symlink(eq_file, joinpath(tmpdir, eq_name))
+    catch
+        cp(eq_file, joinpath(tmpdir, eq_name))
+    end
+    open(joinpath(tmpdir, "gpec.toml"), "w") do io
+        println(io, """
+        # Auto-generated by benchmark_a10_kinetic_dcon.jl — do not edit.
+
+        [Equilibrium]
+        eq_filename = "$eq_name"
+        eq_type = "efit"
+        jac_type = "hamada"
+        power_bp = 0
+        power_b = 0
+        power_r = 0
+        grid_type = "ldp"
+        psilow = 1e-3
+        psihigh = 0.99
+        mpsi = 16
+        mtheta = 256
+        newq0 = 0
+        etol = 1e-7
+
+        [Wall]
+        shape = "nowall"
+
+        [ForceFreeStates]
+        bal_flag = false
+        mat_flag = true
+        ode_flag = true
+        vac_flag = true
+        mer_flag = true
+        force_termination = true
+
+        set_psilim_via_dmlim = false
+        psiedge = 1.0
+        qlow = 1.02
+        qhigh = 1e3
+        sing_start = 0
+
+        nn_low = 1
+        nn_high = 1
+        delta_mlow = 6
+        delta_mhigh = 6
+        delta_mband = 0
+        mthvac = 512
+        thmax0 = 1
+
+        kinetic_source = "calculated"
+        kinetic_factor = 1.0
+        eulerlagrange_tolerance = 1e-7
+        singfac_min = 1e-4
+        ucrit = 1e4
+        write_outputs_to_HDF5 = true
+
+        [KineticForces]
+        kinetic_file = "$kin_file"
+        nn = 1
+        nl = 6
+        zi = 1
+        mi = 2
+        zimp = 6
+        mimp = 12
+        electron = false
+        nutype = "harmonic"
+        f0type = "maxwellian"
+        atol_xlmda = 1e-9
+        rtol_xlmda = 1e-5
+        atol_psi = 1e-9
+        rtol_psi = 1e-4
+        """)
+    end
+    return tmpdir
+end
+
+_p(args...) = (println(stderr, args...); flush(stderr))
+_pf(fmt, args...) = (print(stderr, Printf.format(Printf.Format(fmt), args...)); flush(stderr))
+
+"""
+    eigenvector_overlap(xi_J, mlow_J, xi_F, mlow_F, mhigh_F) → Float64
+
+Align Julia and Fortran total-energy eigenvectors on the intersection of
+their m-ranges, then return `|⟨ξ_J|ξ_F⟩| / (‖ξ_J‖·‖ξ_F‖)`. Vectors are
+phase-insensitive and normalized independently — a value of 1.0 means
+identical shape.
+"""
+function eigenvector_overlap(xi_J::AbstractVector, mlow_J::Int,
+                             xi_F::AbstractVector, mlow_F::Int, mhigh_F::Int)
+    mhigh_J = mlow_J + length(xi_J) - 1
+    m_lo = max(mlow_J, mlow_F)
+    m_hi = min(mhigh_J, mhigh_F)
+    m_lo <= m_hi || return 0.0
+    iJ = (m_lo - mlow_J + 1):(m_hi - mlow_J + 1)
+    iF = (m_lo - mlow_F + 1):(m_hi - mlow_F + 1)
+    vJ = @view xi_J[iJ]
+    vF = @view xi_F[iF]
+    denom = norm(vJ) * norm(vF)
+    denom == 0 && return 0.0
+    return abs(dot(vJ, vF)) / denom
+end
+
+function run_benchmark(fortran_kin_dir::String=DEFAULT_FORTRAN_KIN_DIR)
+    _p("=" ^ 70)
+    _p("  a10 Kinetic-DCON Benchmark: Julia KF→FFS vs Fortran DCON")
+    _p("  Fortran example: $fortran_kin_dir")
+    _p("=" ^ 70)
+
+    inputs = discover_inputs(fortran_kin_dir)
+    ref    = read_dcon_reference(inputs.dcon_nc)
+    _pf("  EFIT g-file: %s\n", basename(inputs.eq_file))
+    _pf("  Kinetic:     %s\n", basename(inputs.kin_file))
+    _pf("  Reference W_t[1]: Re = %.6f  Im = %.6f   (mpert=%d, m∈[%d, %d])\n",
+        ref.W_t_re, ref.W_t_im, ref.mpert, ref.mlow, ref.mhigh)
+
+    tomldir = build_benchmark_tomldir(inputs.eq_file, inputs.kin_file)
+
+    _p("\n--- Running GPE.main with kinetic_source=\"calculated\" ---")
+    t0 = time()
+    GPE.main([tomldir])
+    wall = time() - t0
+    _pf("  GPE.main completed in %.1f s\n", wall)
+
+    h5path = joinpath(tomldir, "gpec.h5")
+    isfile(h5path) || error("Expected Julia output not found: $h5path")
+
+    et_J, xi_J, mlow_J = h5open(h5path, "r") do h5
+        (read(h5["vacuum/et"]),
+         read(h5["vacuum/et_eigenvector"]),
+         read(h5["info/mlow"]))
+    end
+    isempty(et_J)  && error("vacuum/et is empty in $h5path")
+    isempty(xi_J)  && error("vacuum/et_eigenvector is empty in $h5path")
+    _pf("  HDF5 shapes: et=%s, et_eigenvector=%s, mlow=%s\n",
+        string(size(et_J)), string(size(xi_J)), string(mlow_J))
+
+    julia_re = real(et_J[1])
+    julia_im = imag(et_J[1])
+    re_err   = abs(julia_re - ref.W_t_re) / abs(ref.W_t_re) * 100
+    im_err   = abs(julia_im - ref.W_t_im) / abs(ref.W_t_im) * 100
+    overlap  = eigenvector_overlap(xi_J[:, 1], Int(mlow_J),
+                                   ref.xi_F, ref.mlow, ref.mhigh)
+
+    _p("\n" * "=" ^ 70)
+    _p("  Results Comparison (kinetic-DCON least-stable eigenvalue + eigenvector)")
+    _p("=" ^ 70)
+    _pf("  Re(et[1]):  Julia = %12.6f  Fortran = %12.6f  err = %6.2f%%\n",
+        julia_re, ref.W_t_re, re_err)
+    _pf("  Im(et[1]):  Julia = %12.6f  Fortran = %12.6f  err = %6.2f%%\n",
+        julia_im, ref.W_t_im, im_err)
+    _pf("  |⟨ξ_J|ξ_F⟩|: %.4f   (mpert_J=%d, mpert_F=%d)\n",
+        overlap, size(xi_J, 1), ref.mpert)
+    _p("\n  Plasma/vacuum partition (Fortran least-stable, separate diag.):")
+    _pf("    W_p[1] = %.6f + %.6fi\n", ref.W_p_re, ref.W_p_im)
+    _pf("    W_v[1] = %.6f + %.6fi\n", ref.W_v_re, ref.W_v_im)
+
+    re_pass      = re_err   <= 10.0
+    im_pass      = im_err   <= 30.0
+    overlap_pass = overlap  >= 0.90
+
+    _p("\n" * "=" ^ 70)
+    if re_pass && im_pass && overlap_pass
+        _p("  PASS — Re<10%, Im<30%, overlap>0.90")
+    else
+        tags = String[]
+        re_pass      || push!(tags, @sprintf("Re err %.1f%% > 10%%",       re_err))
+        im_pass      || push!(tags, @sprintf("Im err %.1f%% > 30%%",       im_err))
+        overlap_pass || push!(tags, @sprintf("overlap %.3f < 0.90",        overlap))
+        _p("  FAIL — " * join(tags, "; "))
+    end
+    _p("=" ^ 70)
+    _pf("  Wall time: %.1f s\n", wall)
+
+    # cond(F̄) vs ψ diagnostic.
+    png_path = joinpath(@__DIR__, "a10_kinetic_cond_fbar.png")
+    plot_cond_fbar_scan(h5path, png_path; title="a10 kinetic-DCON cond(F̄) vs ψ")
+
+    return (; re_err, im_err, overlap, wall, julia_re, julia_im,
+              fortran_re=ref.W_t_re, fortran_im=ref.W_t_im)
+end
+
+if abspath(PROGRAM_FILE) == @__FILE__
+    fortran_dir = length(ARGS) >= 1 ? ARGS[1] :
+                  get(ENV, "GPEC_FORTRAN_A10_DCON", DEFAULT_FORTRAN_KIN_DIR)
+    run_benchmark(fortran_dir)
+end