API Reference

Per-radius marginal TGLF-EP QL mode data (one scan point).

Options{T<:Real}

Scan-control / run settings for a TJLFEP calculation — the Julia counterpart of the Fortran TGLFEP interface module populated from an input.TGLFEP file.

Holds the process/threshold selectors (PROCESS_IN, THRESHOLD_FLAG, SCAN_METHOD, rejection-flag toggles), basis/mode sizes (N_BASIS, NMODES, NTOROIDAL), the radial scan definition (SCAN_N, IRS, IR_EXP), the EP/ion species descriptors (M_I/Q_I/M_EP/Q_EP, IS_EP), the width controls (WIDTH_IN_FLAG, WIDTH_IN/WIDTH_MIN/WIDTH_MAX), and the per-(kyhat,width) keep/reject result masks (LKEEP, LTEARING, the pinch flags, ...) plus scan outputs (FACTOR_OUT, WIDTH, F_REAL).

Fields default to missing/NaN and are filled by the readers (readTGLFEP, TJLFEP_generate_input); the inner constructor Options{T}(nscan_in, widthin, nn, nr, jtscale_max, nmodes) allocates the vector/matrix fields to the right sizes.

TJLFEP_TM(inputsEP::Options, inputsPR::profile; mode_in, use_gpu=false)

TGLFEP spectrum diagnostic (PROCESS_IN=3) port. Computes the gyro-Bohm normalized growth-rate/frequency spectra γ(ky), ω(ky) by running the full TJLF transport model over a fixed ky grid (nky=30, ky=0.15, KYGRID_MODEL=0, NBASIS=32), matching the Fortran TGLFEP_TM / write_eigenvalue_spectrum.

mode_in selects the drive (via TJLF_map's override): 1 = background plasma + EPs, 2 = EP-only, 4 = ITG/TEM with TAE/EPM filtered. Pass use_gpu=true to run the eigensolves on CUDA.

Returns (ky, gamma, freq, (filename, buffer)), where gamma/freq are nky × nmodes matrices (row i = ky[i]) and buffer reproduces out.eigenvalue_m<mode><suffix>.

TJLFEP_generate_input()

Generate the TJLFEP profile input files from an input.TJLFEP.generate source in the current directory (a port of the Fortran TGLFEP_generate_input utility).

Reads the sign/geometry/species header and per-radius profile columns and writes the corresponding TJLFEP input(s). Errors out if input.TJLFEP.generate is not found.

TJLFEP_ky_widthscan(inputsEP::Options, inputsPR::profile; use_gpu=false, inner=:threads, team=nothing)

TGLFEP ky width-scan port, used by the spectrum diagnostic (TJLFEP_TM, PROCESS_IN=3) when WIDTH_IN_FLAG=false. Scans the single-ky growth-rate/frequency over a width grid (WIDTH_MIN : 0.01 : WIDTH_MAX) at the EP-only drive and picks the width that maximizes the kept-mode growth rate.

inner/team select the parallel backend (:threads, or :mps_team with GPU ids); use_gpu=true runs eigensolves on CUDA.

Returns the chosen width, the ky used (constant across the scan → kymark), and the out.ky_widthscan_m<mode> buffer.

TJLF_map(inputsEP::Options, inputsPR::profile; mode_in_override=nothing, ky_model_override=nothing, nky=12) -> TJLF.InputTJLF

Build a TJLF.InputTJLF for a single radial point from TJLFEP's Options (scan control) and profile (plasma profiles), mirroring the Fortran TGLFEP → TGLF variable mapping.

mode_in_override / ky_model_override force the drive selector and ky-grid model (used by the spectrum diagnostic), and nky sets the ky-grid size.

build_alpha_ql_modes(marginals, rho_scan, rho_grid, dndr_crit; km_max=5) -> Vector

Map per-scan MarginalQLData to a radial vector of named tuples for ALPHA.QLModeInput.

compute_alpha_crit_profiles(SFmin, Options, profile, ni, Ti, dlnnidr, dlntidr)

Recompute alpha_dndr_crit / alpha_dpdr_crit profiles from converged SFmin (same logic as run_tjlfep_file).

critical_factor_ad_f1seed(ep0, prof; kwargs...)

Cheap "front-end" critical-factor solver: an f1-seed grid over (kyhat, width) followed by an AD (:ad) descent on the interior basins and an early-stop few-confirm pass. Floor-pinned basins (onset at the scan floor) are caught via the AE-unstable-window lower-edge f1, while interior basins are slid off-node by the descent. This is the fast tier feeding the :ad :locate/robust_ad escalation ladder.

Key keywords: gamma_thresh (growth-rate threshold; defaults per case), scan_lo/scan_hi (factor bracket; default FACTOR_IN and FACTOR_IN/512), nseed_ky/nseed_w/n_eig_seed (seed-grid sizes), k_descend (interior basins to descend), ky_lo=0.25 (kyhat floor), and inner/team/use_gpu for the parallel/GPU backend.

Returns the critical factor together with trust diagnostics (e.g. feasible_frac, cheap_gap) used by the escalation gate.

critical_factor_confirm(inputsEP, inputsPR; kwargs...) -> NamedTuple

Cheap-search + few-confirm critical EP scale factor. Scans the (kyhat, width) grid on the CHEAP eigenvalue-only AE-band onset (IFLUX=false), then applies the expensive faithful keep filters (IFLUX=true, the tearing/pinch/QL-ratio/θ² block) at ONLY the few candidate nodes that could hold the global minimum.

This is exact, not heuristic: the faithful keep onset at any (ky,w) is ≥ its AE-band onset (keep ⊆ AE-unstable, since keep requires γ>γ* AND in-band AND surviving the secondary filters). So candidates are confirmed in increasing cheap-onset order and the loop stops as soon as the next candidate's cheap onset exceeds the best faithful onset found — no remaining node can beat it. The result therefore EQUALS the faithful (ky,w)-grid minimum (the critical_factor_robust refine=0 answer), but pays IFLUX=true on only ~1–few nodes instead of all nkyhat·nefwid. The saving is large when the secondary filters rarely bind (e.g. DIII-D, where the binding is AE-band growth and the pinch flags are off); when they bind often (rotational ITER) more nodes are confirmed, degrading gracefully to the full faithful grid.

The (ky,w) grid mirrors kwscale_scan's first round (kyhat[i]=(1/nkyhat)·i, width[i]=WIDTH_MIN+(WIDTH_MAX-WIDTH_MIN)/(nefwid-1)·(i-1)). Keywords: nkyhat=4, nefwid=8, gamma_thresh=nothing (case γ*), scan_lo=nothing(→scan_hi/(nfactor·4^(k_max-1)), the grid floor; scan_hi/512 at kmax=4), `scanhi=nothing(→FACTORIN),inner=:threads,team=nothing,usegpu=false,verbose=false`.

Returns (; sfmin, kyhat, width, binding, status, cheap_min, n_confirm, npts, total_evals_full, total_evals_eig, results), where status ∈ (:ok, :no_onset, :cap), cheap_min is the lowest cheap AE-band onset over the grid, and n_confirm is the number of nodes that actually got the IFLUX=true confirm.

critical_factor_grid(inputsEP, inputsPR; kwargs...) -> NamedTuple

AD analogue of the inner kwscale_scan: the critical EP scale factor is the minimum over the (kyhat, width) grid of the per-point marginal factor. For each grid point this calls marginal_factor (AD-Newton) instead of the brute-force 8-factor × 4-round bracket, so it reaches sfmin in far fewer TGLF eigensolves.

The (kyhat, width) grid mirrors kwscale_scan's first round:

• kyhat[i] = (1/nkyhat)·i, i=1..nkyhat (kyhat∈(0,1])
• width[i] = WIDTH_MIN + (WIDTH_MAX-WIDTH_MIN)/(nefwid-1)·(i-1)

Keywords: nkyhat=4, nefwid=8, gamma_thresh=1e-7, scan_lo=1e-3, scan_hi=nothing (defaults to the incoming FACTOR_IN), threaded=true, use_gpu=false.

Returns (; sfmin, kyhat, width, gamma, total_evals, results) where results is the per-grid-point table.

critical_factor_multistart(inputsEP, inputsPR; local_algo=:LN_BOBYQA, kwargs...) -> NamedTuple

Derivative-free (kyhat, width) critical-factor search: cheap AE-onset multistart + local derivative-free descents on the cheap onset surface + faithful early-stop confirm. The bumpy onset surface f1(ky,w) is multimodal, so global-ness comes from the multistart (_cheap_onset_seeds); the local driver only refines within a basin:

1. Cheap eigenvalue-only (IFLUX=false) AE-onset rank over a (log+linear width, down to w_lo) × kyhat grid → feasible seeds sorted by onset.
2. For the top k_descend well-separated non-pinned seeds, run a local derivative-free minimization of the cheap onset in (ky, log w) coordinates (log-width matches the band's natural scaling). local_algo selects the driver: :LN_BOBYQA (default) — Powell quadratic trust-region model; most eval-efficient. :LN_SBPLX — Rowan subplex; noise-robust on the bumpiest surfaces. :LN_NELDERMEAD — plain simplex. Every point the optimizer visits is recorded, and each cheap eval is ONE batched-parallel _ae_unstable_window (~10-20× cheaper than a marginal_factor_df f★ solve).
3. Faithful-confirm (IFLUX=true) the pooled candidates — ALL feasible cheap seeds (floor guards) plus all unstable descent samples — in ascending cheap order with an early-stop bound (_confirm_candidates) → the critical factor.

Keywords: local_algo=:LN_BOBYQA, gamma_thresh=nothing, scan_lo=nothing(→shi/(nfactor·4^(k_max-1)); shi/512 at kmax=4), `scanhi=nothing(→FACTORIN),kylo=0.05,wlo=0.05,nseedky=6,nseedw=10,neigseed=12,neig=16,kdescend=3,localevals=15,maxconfirm=8,inner=:threads,team=nothing,usegpu=false,verbose=false. The faithful confirm is the wall-time driver:max_confirm` hard-caps it (the early-stop prune plus basin dedup usually stop well before the cap).

Returns (; sfmin, kyhat, width, binding, status, converged, faithful, n_confirm, n_samples, total_evals_full, total_evals_eig), status ∈ (:ok,:no_onset,:cap).

critical_factor_nlopt(inputsEP, inputsPR; kwargs...) -> NamedTuple

Derivative-free (kyhat, width) critical-factor search using NLopt global derivative-free optimization on the cheap AE-band onset surface, then a faithful early-stop confirm. This is the core promotion of the offline critical_factor_direct prototype, generalized to a selectable algorithm and an optional local polish:

1. NLopt global search (algo, default :GN_DIRECT_L; also :GN_CRS2_LM, :GN_AGS, ...) over (ky,w) minimizing the cheap AE onset f1 (_ae_unstable_window, IFLUX=false; stable points penalized to FACTOR_IN), across the narrow-width-extended box ky∈[ky_lo,1], w∈[w_lo,WIDTH_MAX].
2. Local polish (local_algo, default :LN_BOBYQA; also :LN_SBPLX/:LN_COBYLA, nothing disables) seeded at the global best — DIRECT-L brackets the winning basin but converges slowly inside it; BOBYQA's quadratic model closes the last stretch in a handful of evals.
3. Faithful-confirm (IFLUX=true) the sampled points in ascending cheap order with early stop.

Keywords: algo=:GN_DIRECT_L, local_algo=:LN_BOBYQA, max_evals=40, local_evals=20, gamma_thresh=nothing, scan_lo=nothing, scan_hi=nothing, ky_lo=0.05, w_lo=0.05, n_eig=16, inner=:threads, team=nothing, use_gpu=false, verbose=false.

Returns the same contract as critical_factor_multistart.

critical_factor_optimize(inputsEP, inputsPR; kwargs...) -> NamedTuple

AD analogue of the inner kwscale_scan: find the critical EP scale factor sfmin = min over (kyhat,width) of the marginal factor f★(kyhat,width) by projected-gradient descent with implicit-function-theorem gradients.

Each evaluated (ky,w) solves the AE-band growth root f★ by AD-Newton (marginal_factor) and gets the exact ∂f★/∂(ky,w) from one extra gamma_grad pass (the implicit-function-theorem ratio -(∂γ/∂x)/(∂γ/∂f)). A coarse nseed_ky × nseed_w seed grid (or a caller-supplied seed) locates a feasible basin; descent then polishes it with backtracking line search, projecting (ky,w) to ky∈ky_bounds, w∈[WIDTH_MIN,WIDTH_MAX].

NOTE on cost: this yields the continuous optimum (more precise than any fixed grid), but it is NOT necessarily cheaper than a (ky,w) grid scan — the AE-band onset surface f★(ky,w) is bumpy/multimodal (the same narrow-bump physics the keep window shows in factor), so descent can need comparably many TGLF solves and warm seeds can land in a worse basin. The large AD win is per-evaluation (eigenvalue-only IFLUX=false Newton vs the production IFLUX=true factor grid), realized in marginal_factor/marginal_factor_faithful, not in collapsing the (ky,w) grid.

Keywords: gamma_thresh=nothing (case γ*), nseed_ky=4, nseed_w=4, scan_lo=1e-3, scan_hi=nothing(→FACTORIN), nscan=8, `kybounds=(1e-3,1.0),wbounds=nothing(→struct),maxiter=25,gtol=1e-4,xtol=1e-5,step0=0.25,nbacktrack=12,faithfulconfirm=false,use_gpu=false,verbose=false`.

extend_width=false (opt-in; on by default in the solver=:ad scan path): the canonical descent above only searches w∈[WIDTH_MIN,WIDTH_MAX] (w≥1), so it structurally misses the narrow EP-driven AE modes (w≪1) that set the true critical factor at near-marginal radii (the canonical w≥1 box biases sfmin high there — up to ~25× at the steep profile edge). When extend_width=true the solver additionally runs _locate_extended — a cheap log-spaced (ky,w) AE-onset rank (down to w_lo=0.05, ky_lo=0.05) → :ad descent → faithful early-stop confirm — and folds any LOWER faithful onset into the result. This is the SAME width-extension machinery :robust_ad uses, but seeded on the cheap :ad descent instead of a full faithful w≥1 grid-zoom, so it is much cheaper than :robust_ad (no faithful grid) while recovering most of the narrow-width accuracy gain over the Fortran/grid w≥1 baseline. The w≥1 descent is retained as a floor (AD descent can under-resolve the dense core). Extension knobs: w_lo=0.05, ky_lo=0.05, nseed_ky_ext=6, nseed_w_ext=10, n_eig_seed=12, k_descend=6. It does NOT alter the default (extend_width=false) behavior, so _locate_extended's own inner :ad descents (which pass it unset) remain w≥1.

extend_mode (default :locate) selects HOW the extension searches:

• :locate — the _locate_extended log-grid + k_descend-way multistart above (production; matches :robust_ad bitwise under faithful_confirm).
• :wide — a much cheaper single-pass scheme: widen the main descent box down to w_lo (and clamp the ky floor up to ky_lo to match :locate's domain), seed it with a denser log-spaced w grid, descend the top wide_kdesc (default 2) well-separated seeds, drop any floor-pinned (w≈w_lo/ky≈ky_lo) corner optima when an interior optimum exists, and faithful-confirm the best. ~7× cheaper than :locate, recovers most of the narrow-width gain over the w≥1 baseline, and stays conservative (confirmed keep onset ≥ truth) — it does NOT generally match :robust_ad bitwise. Good for fast NN-database generation. With faithful_confirm=false :wide reports the (unfiltered) AE-band onset and can under-predict, like the pure :locate path.

extend_width interacts with faithful_confirm:

• faithful_confirm=true (the solver=:ad scan default): the located narrow candidates are faithful-confirmed (IFLUX=true) with an early-stop bound, so sfmin is the faithful keep onset (min over the w≥1 descent and the extension) and matches :robust_ad bitwise at the located mode.
• faithful_confirm=false ("pure AD" extended): NO faithful confirm anywhere — the solver folds in the lowest cheap AE-band/descended onset over the located candidates and reports that. Much cheaper (zero IFLUX=true evals), but the AE-band onset is ≤ the faithful keep onset, so it does not match :robust_ad bitwise (and can dip below it where keep filters would reject the marginal mode).

Returns (; sfmin, kyhat, width, f_seedmin, ky_seed, w_seed, iters, evals, converged, faithful, extended, n_ext_confirm) where faithful is the marginal_factor_faithful result at the winning optimum (or nothing), extended echoes extend_width, and n_ext_confirm counts the narrow-width faithful confirms (0 when faithful_confirm=false). With a faithful confirm, sfmin is the faithful onset of the winner (min over the w≥1 descent and the narrow extension); otherwise the AE-band descent onset.

critical_factor_profile(inputsEP, inputsPR; kwargs...) -> NamedTuple

AD analogue of the radial driver loop (tjlfep_driver / mainsub over the SCAN_N radii): the critical EP scale factor sfmin(IR) versus radius, computed by critical_factor_optimize at each radius with cross-radius continuation — the optimum (kyhat,width) of one radius warm-starts the next. (Continuation skips the per-radius seed grid, but on the bumpy f★(ky,w) surface it does not reliably reduce total solves and can land in a worse local basin; see the cost note on critical_factor_optimize. It is most useful when the optimum moves smoothly with radius.)

Keywords: radii=nothing (→ Int.(IR_EXP[1:SCAN_N])), scan_lo=1e-3, scan_hi=10.0, faithful=false (confirm each optimum with marginal_factor_faithful), use_gpu=false, verbose=false; remaining kwargs forwarded to critical_factor_optimize.

Returns (; radii, sfmin, kyhat, width, evals, converged, binding, total_evals). When faithful=true, sfmin holds the faithful (all-filter) onset and binding the per-radius binding condition; otherwise sfmin is the AE-band onset.

critical_factor_robust(inputsEP, inputsPR; kwargs...) -> NamedTuple

Robust, grid-faithful critical EP scale factor: the global minimum over a (kyhat, width) grid of the FAITHFUL (all-filter) marginal factor, mirroring exactly what the Fortran/grid kwscale_scan reduces to (sfmin = min over (kyhat,width) of the kept-window lower edge). Unlike critical_factor_optimize, which gradient-descends the bumpy/multimodal AE-band onset surface f★(ky,w) and can land in a shallower local basin (or on a mode the keep filters later reject), this evaluates marginal_factor_faithful at every grid point and takes the global min — so it reproduces the grid sfmin to within the continuous-vs-discrete-factor difference, at the cost of scanning the full (ky,w) grid (no (ky,w) collapse). The inner factor root is still the cheap AD-Newton + faithful-bisect (IFLUX=false hull localization, IFLUX=true only inside the AE band), so each grid point is far cheaper than the production 8-factor × 4-round brute force.

The coarse (kyhat, width) grid mirrors kwscale_scan's first round:

• kyhat[i] = (1/nkyhat)·i, i=1..nkyhat
• width[i] = WIDTH_MIN + (WIDTH_MAX-WIDTH_MIN)/(nefwid-1)·(i-1)

refine_rounds then re-grids a shrinking (ky,w) box around the running best, mirroring the Fortran window narrowing so binding points that lie between coarse nodes (e.g. the plasma edge) are resolved rather than overestimated. With adaptive=true (default), refine_rounds is a per-radius budget: a round runs only when an a-priori trigger fires — the best is on a box edge (minimum unbracketed), the feasible region is sparse (≤2 or < feasible_frac of nodes carried a kept mode), or the best onset is near the cap (≥ cap_frac·scan_hi) — and stops early once a round's relative improvement drops below improve_tol and the best is interior. So radii whose coarse min is already interior/dense/low spend no refinement, while edge/sparse radii get up to refine_rounds zoom rounds. With adaptive=false, exactly refine_rounds rounds always run.

outer selects how the off-node minimum is resolved once the coarse grid has bracketed the basin:

• :grid_zoom (default) — re-grid a shrinking (ky,w) box (the refine_rounds loop above);
• :hybrid — when the zoom trigger is QUIET (best interior, feasibility dense, onset below cap), replace the grid-zoom round with a LOCAL AD-IFT polish (critical_factor_optimize seeded at the coarse best: projected-gradient on the cheap IFLUX=false AE-onset surface using exact ∂f★/∂(ky,w)) followed by a FAITHFUL confirm of the optimum. The coarse grid still does the (multimodality-safe) global bracketing; AD only does the local descent it is good at. When the trigger FIRES (untrustworthy bracket), :hybrid falls back to the grid-zoom loop. This aims to recover the off-node accuracy of a refine round at a handful of solves instead of a full grid.

extend_width (default true) extends the search BELOW the canonical box to capture the radially narrow EP-driven AE modes (width ≪ WIDTH_MIN) that w∈[WIDTH_MIN,WIDTH_MAX] structurally misses. After the canonical grid-zoom (which is kept as a robust FLOOR), it runs _locate_extended over a log-spaced width box down to w_lo — cheap IFLUX=false AE-onset ranking → :ad (AD-IFT) descent of the top seeds → FAITHFUL confirm — and folds any lower faithful onset into the global best. So sfmin = min(canonical w≥1 grid-zoom, extended narrow-width locate). This is the width component of the grid→truth gap and, on the DIII-D edge, drops sfmin by 2–11× vs the w≥1-only value. Set extend_width=false to recover the pure faithful w≥1 minimum (≈ the :grid reference). The nbasis-converged tier is critical_factor_truth, which is this solver plus a nbasis ladder at the located optimum.

confirm_grid (default false) is an OPT-IN cheap-rank→few-confirm scheme for the w≥1 grid passes: rank every (ky,w) node on the cheap IFLUX=false AE-band onset and faithful-confirm only the few that can still hold the minimum (the critical_factor_confirm logic, shared via _rank_confirm). It is EXACT for a fixed node set (verified bitwise with adaptive=false) and cuts IFLUX=true evals, but is kept OFF by default because on the GPU/MPS production route it is measurably SLOWER (the early-stop confirm is serial — one GPU — while the brute path fans all nodes across the MPS team; the keep filter is not the per-node bottleneck the eigensolve is) and, with adaptive=true, its cheap (over-counted) feasibility count shifts the sparse zoom trigger so sfmin can differ from the brute path at near-degenerate radii. Prefer it for CPU/serial runs or with adaptive=false.

Keywords: nkyhat=4, nefwid=8, refine_rounds=1, adaptive=true, confirm_grid=false, outer=:grid_zoom, feasible_frac=0.35, cap_frac=0.5, improve_tol=0.02, extend_width=true, ky_lo=0.05, w_lo=0.05, nseed_ky=6, nseed_w=10, n_eig_seed=12, k_descend=6, gamma_thresh=nothing (case γ*), scan_lo=nothing(→scan_hi/(nfactor·4^(k_max-1)), the grid floor; scan_hi/512 at kmax=4), `scanhi=nothing(→FACTORIN),inner=:threads,team=nothing,usegpu=false,verbose=false. The(ky,w)` points are parallelized over the team/threads; each point's inner faithful sweep runs serially to avoid nested oversubscription.

Returns (; sfmin, sfmin_w1, kyhat, width, binding, status, scan_lo, scan_hi, refine_done, polished, extended, n_ext_confirm, n_feasible_coarse, npts_coarse, total_evals_full, total_evals_eig, npts, results), where status ∈ (:ok, :no_onset, :cap), sfmin_w1 is the canonical w≥1 grid-zoom minimum before the width extension (the well-converged Fortran-faithful threshold, used by critical_factor_truth as a floor), refine_done is the number of grid-zoom rounds performed, polished is true when the :hybrid AD-IFT polish ran, extended reflects extend_width, and n_ext_confirm is the number of faithful confirms spent on the width extension.

critical_factor_triggered(inputsEP, inputsPR; w_floor_tol=0.05, gap_thresh=1.5, feas_thresh=0.25,
                          ky_lo=0.25, truth_kwargs=NamedTuple(), inner=:threads, team=nothing,
                          use_gpu=false) -> NamedTuple

Production policy: run the fast canonical critical_factor_ad_f1seed pass on w∈[WIDTH_MIN,WIDTH_MAX], and escalate to critical_factor_truth only when the canonical optimum pins at the width floor (a narrow-width bowl likely lies below the box) or the trust gate fires (sparse feasible fraction, cheap_gap keep-filter divergence, no_onset, cap). When escalated, returns the lower (min(canonical, truth)) faithful onset. Returns (; sfmin, kyhat, width, binding, status, escalated, reasons, source, base_sfmin, truth_sfmin, ...).

Note: because narrower width is systematically more unstable, the canonical optimum pins at WIDTH_MIN at most radii, so the :width_floor trigger is broad. Use truth_kwargs to forward solver knobs (e.g. nseed_w, nb_steps) to the escalated search.

critical_factor_truth(inputsEP, inputsPR; ky_lo=0.05, w_lo=0.05, nseed_ky=6, nseed_w=10,
                      nb_work=32, nb_steps=[32,40,48,56], inner=:threads, team=nothing,
                      use_gpu=false) -> NamedTuple

The nbasis-converged accuracy tier. This is critical_factor_robust(extend_width=true) plus a nbasis ladder at the located optimum — the third rung of the grid → robust_ad → truth cost/ fidelity ladder. It first calls critical_factor_robust (width-extended: min of the canonical w≥1 grid-zoom and the extended narrow-width locate) to find the production (ky,w) and its faithful sfmin at nb_work, then re-evaluates the FAITHFUL marginal factor at that FIXED optimum over nb_steps (geometric extrapolation, capped ≤56; SingularException at high (nb,narrow-w) → NaN, not fatal). The nbasis correction is adverse (raises the threshold) at the still-converging outer radii — exactly the +nbasis component of the grid→truth gap (the width component, the 2–11× reduction, is already captured by robust_ad). The reported sfmin is the nbasis-converged value floored by robust_ad's canonical w≥1 minimum (sfmin_w1): a razor-thin narrow win at nb_work can flip above the w≥1 threshold once converged, and the true critical factor is the lowest threshold over ALL modes (the w≥1 mode is already converged at nb_work), so sfmin = min(nbasis-converged narrow, w≥1). On most radii the narrow mode stays well below w≥1 and the floor is a no-op. Returns (; sfmin, sfmin_work, sfmin_conv, sfmin_w1, floored, kyhat, width, binding, status, nb_limit, nb_ratio, nb_converged, nb_table, feasible_frac, n_confirm, total_evals_full, total_evals_eig) where sfmin is the floored production value, sfmin_conv the raw nbasis-converged narrow value before the floor, sfmin_work the robust_ad value at nb_work, and floored whether the w≥1 floor won.

diff_star_from_D_W(ep_flux, factors, inputsEP, inputsPR) -> Vector

Fortran: D_W_k = (dFlux_W_k/dn_sf) / (-d n_sd/dr) with flux linear in n_sf.

expro_bound_deriv(f, r) -> Vector{Float64}

Match GACODE expro_util.f90::bound_deriv (Lagrange derivative on uniform radial grid).

Build EXPRO dict for save_EXPRO from a profile (3-species Fortran layout).

dlnnidr = -d(ln n)/dr and dlntidr = -d(ln T)/dr as in expro_compute_derived.

Map GACODE is_EP (ions only) to EXPRO species index when species 1 is electron.

expro_vectors_from_gacode(prof, gacode_file, is_EP_gacode)

EP α-postprocessing vectors from input.gacode only (no input.EXPRO file). is_EP_gacode is the GACODE ion index from input.TGLFEP (IS_EP).

extract_marginal_ql(inputsEP, inputsPR, growthrate, dep_scan, factor; ...) -> MarginalQLData

gamma_star from the stability scan; diff_star from the Fortran flux scan only when use_flux_scan=true.

finalize_gacode_scan(gacode_file, tglfep_file, tasks_dir; printout)

Merge per-radius task_*.jls files and write alpha_*_crit.input (and optional out.TGLFEP header).

gamma_dgamma_dfactor(inputsEP, inputsPR; use_gpu=false)

Single-point AD building block for the EP critical-gradient scan.

Evaluates the TGLF-EP dispersion at one (KYHAT_IN, WIDTH_IN, FACTOR_IN, IR) operating point — exactly as the inner kwscale_scan combo does — and returns both the per-mode growth rate γ and its exact derivative dγ/d(FACTOR_IN), computed in a single forward-mode AD pass by seeding a ForwardDiff.Dual on the EP scale factor.

The derivative is taken through TJLF_map (by promoting Options/profile to a Dual element type and seeding FACTOR_IN), so every channel by which the factor reaches the TGLF input is captured automatically — not only the EP density gradient RLNS[is] / density AS, but also the pressure-gradient (P_PRIME_LOC, MHD-α) stabilization term, which for PPRIME_METHOD ∈ {2,3} varies with the scaled EP gradient via the sum*/sum1 ratio. Hand-replicating the factor→input map is brittle for exactly this reason.

Returns a NamedTuple:

• gamma :: Vector{Float64} per-mode growth rate (length NMODES)
• dgamma_dfactor :: Vector{Float64} per-mode dγ/d(FACTOR_IN)
• freq :: Vector{Float64} per-mode real frequency
• dfreq_dfactor :: Vector{Float64} per-mode d(freq)/d(FACTOR_IN)
• factor :: Float64 the (clamped) factor the result is evaluated at

inputsEP.IR, FACTOR_IN, KYHAT_IN, and WIDTH_IN must be set on entry (as the scan sets them). inputsEP/inputsPR are not mutated.

gamma_grad(inputsEP, inputsPR; vars=(:FACTOR_IN,:KYHAT_IN,:WIDTH_IN), use_gpu=false)

Multi-variable generalization of gamma_dgamma_dfactor: one forward-mode AD pass returning the per-mode growth rate γ and real frequency together with their exact partials w.r.t. each operating-point variable in vars — any subset of (:FACTOR_IN, :KYHAT_IN, :WIDTH_IN). All partials are obtained in a single N-partial ForwardDiff.Dual solve (one eigensolve), which is the building block for the (kyhat, width) optimizer (Phase 3) and the freq/keep-condition root finds (Phase 2).

Derivatives are taken through TJLF_map (by promoting Options/profile to the Dual element type and seeding the chosen EP fields), so every channel by which a variable reaches the TGLF input is captured automatically. IFLUX=false, so only the eigenvalue block runs.

Returns a NamedTuple:

• gamma :: Vector{Float64} per-mode growth rate (length NMODES)
• freq :: Vector{Float64} per-mode real frequency
• dgamma :: Matrix{Float64} NMODES × length(vars), ∂γ_m/∂vars[k]
• dfreq :: Matrix{Float64} NMODES × length(vars), ∂freq_m/∂vars[k]
• vars :: NTuple the variable order matching the columns

inputsEP.IR, FACTOR_IN, KYHAT_IN, WIDTH_IN must be set on entry (as the scan sets them). inputsEP/inputsPR are not mutated.

gamma_input_sensitivities(inputsEP, inputsPR; knobs=nothing, use_gpu=false) -> NamedTuple

AD-exact sensitivities of the per-mode growth rate γ (and real frequency) to the local TGLF plasma inputs, at a fixed operating point (IR, KYHAT_IN, WIDTH_IN, FACTOR_IN). This is the genuinely differentiable, grid-independent quantity the Fortran TGLF-EP cannot provide — unlike the critical factor sfmin, which sits at a discrete keep-flag transition and is not smoothly differentiable.

The growth rate is differentiated through TGLF's eigensolve in a single forward pass by seeding multi-partial ForwardDiff.Duals on the selected inputs of the InputTJLF produced by TJLF_map (so the eigenvalue rule in TJLF's TJLFForwardDiffExt provides exact ∂λ/∂input). Uses IFLUX=false (eigenvalues only), so it is cheap.

knobs is a vector of (field::Symbol, idx) pairs; idx::Int selects a species entry of a vector field (:RLNS,:RLTS,:TAUS,:AS,:VPAR,:VPAR_SHEAR) and idx=nothing a scalar field (:Q_LOC,:Q_PRIME_LOC,:RMIN_LOC,...). Defaults to per-species a/Ln, a/LT, density, plus Q_LOC.

Returns (; gamma, freq, dgamma, dfreq, logsens, knobs, labels, base) where dgamma/dfreq are NMODES × Nknob Jacobians, base the knob base values, and logsens[m,j] = base[j]·∂γ_m/∂knob_j the dimensionless (logarithmic) sensitivity used to rank drive (>0) vs damping (<0) terms.

ir_exp_from_scan(nr::Int, irs::Int, scan_n::Int) -> Vector{Int}

Compute the scan_n experimental radial indices for a scan: evenly spaced from the start index irs to nr. With scan_n == 1 the single index irs is returned.

make_crit_grad_plots(SorF="neither"; dir="", scale="identity", code="julia")

Plot the critical EP pressure- and density-gradient profiles (alpha_dpdr_crit, alpha_dndr_crit) written by a TJLFEP scan, saving *_dpdr_crit.png / *_dndr_crit.png.

SorF selects the input/output directory and filename prefix ("STRUCT" → structOutputs/, "FILE" → fileOutputs/, anything else → dir/unprefixed). scale sets the y-axis scale and code selects the input file format ("julia" array literal vs. Fortran one-value-per-line).

marginal_factor(inputsEP, inputsPR; kwargs...) -> NamedTuple

AD-accelerated Newton root-find for the marginal EP scale factor at a fixed (IR, KYHAT_IN, WIDTH_IN) operating point: the factor f★ at which the leading growth rate crosses the threshold,

γ_lead(f★) = gamma_thresh .

This is the Julia-native, derivative-based complement to the Fortran 1-D factor scan: instead of bracketing on an 8-point factor grid and taking a single secant step, it runs a safeguarded Newton (Newton step with a bisection fallback) using the exact dγ_lead/d(factor) from gamma_dgamma_dfactor, so it converges quadratically and needs only a handful of TGLF eigensolves.

Keywords:

• gamma_thresh = nothing threshold; nothing ⇒ the case's own GAMMA_THRESH (rotational γ* if ROTATIONAL_SUPPRESSION_FLAG=1, else 1e-7)
• bracket = nothing (f_lo, f_hi) bracketing the sign change of g(f)=γ_lead(f)-gamma_thresh; auto-found if nothing
• scan_lo=1e-3, scan_hi=10.0, nscan=8 geometric samples used to auto-bracket
• ae_band = false track only the AE-band-filtered γ (freq < FREQ_AE_UPPER) instead of the raw max-over-modes γ
• monotone_earlyout = true (raw γlead only; auto-disabled when `aeband) assume γ_lead↑ infand probescanhi/scanlofirst, settling "always stable/unstable" (the high-γ*/ExB case) in ≤2 eigensolves instead ofnscan`
• tol = 1e-6 converged when |g| < tol·(1+|γ_thresh|) or step < xtol
• xtol = 1e-8, maxiter = 40
• use_gpu = false, verbose = false

Returns (; factor, gamma_lead, gamma, freq, iters, evals, converged, bracket).

marginal_factor_df(inputsEP, inputsPR; kwargs...) -> NamedTuple

Derivative-FREE analogue of marginal_factor: the marginal EP scale factor f★ at which the AE-band leading growth rate crosses threshold,

γ_AE(f★) = gamma_thresh ,

found by bracketing + a derivative-free root solver (SimpleNonlinearSolve ITP/Brent/Bisection) instead of the safeguarded Newton. This is the direct replacement for the :grid engine's coarse 8-point factor scan + single linear interpolation (see the "How :grid homes in on the threshold" analysis): it converges to the true threshold at arbitrary precision in a handful of eigensolves and — because it only ever brackets the up-crossing — it cannot overshoot the narrow AE γ_AE(f) bump the way a monotone-assuming interpolation can.

It shares marginal_factor's bracketing preamble (a warm f_start geometric expansion for neighbor-radius continuation, else a log-spaced nscan scan for the first γ_AE(f)-γ* up-crossing), but evaluates only the growth rate (the AD dγ/df is never used by the root solver). Defaults to ae_band=true since the derivative-free solvers track the AE band.

Keywords mirror marginal_factor: gamma_thresh, bracket, f_start, scan_lo=1e-3, scan_hi=10.0, nscan=8, ae_band=true, plus rootsolver=:itp (:itp | :brent | :bisection), tol=1e-6, maxiter=100, use_gpu=false, verbose=false.

Returns the same shape as marginal_factor: (; factor, gamma_lead, gamma, freq, iters, evals, converged, bracket).

marginal_factor_faithful(inputsEP, inputsPR; kwargs...) -> NamedTuple

Ground-truth marginal EP scale factor at a fixed (IR, kyhat, width) point that matches the production keep definition (all secondary IFLUX=true filters), using cheap AD eigensolves to localize the search and full evals only where needed.

Strategy:

1. Cheap eigenvalue-only scan (IFLUX=false) maps the AE-band UNSTABLE hull [f1,f2] where γ_AE(f) > γ*, with f1 Newton-refined. Because the keep window is a subset of this hull (keep ⇒ g>γ*), and γ_AE(f) is typically a narrow bump, restricting the expensive search to [f1,f2] is both complete and immune to stepping over the window.
2. Faithful keep sweep (IFLUX=true) on a fine grid within [f1,f2]: the lowest kept sample is bisected against its unkept neighbor to give the lower edge of the kept window. The condition that fails just below that edge names the binding filter (band-entry/growth for DIII-D; ion/thermal pinch or QL-ratio for the rotational ITER case).

Keywords: kyhat, width (default to the struct's), gamma_thresh=nothing (case γ*), scan_lo=1e-3, scan_hi=10.0, n_eig=24 (eigenvalue scan points), n_fine=28 (faithful samples within the hull), xtol=1e-4, maxbisect=30, use_gpu=false, verbose=false.

Returns (; factor_faithful, factor_fast, window, binding, kept_modes, evals_eig, evals_full, converged, keep). factor_fast is the cheap AD instability onset f1; window=(f1,f2). binding ∈ {:ae_band_growth, :tearing, :i_pinch, :e_pinch, :th_pinch, :ep_pinch, :ql_ratio, :theta_sq, :none}.

populate_tjlfep_profile!(prof::profile, extraEP::Dict, nr::Int, ns::Int)

Populates TJLFEP profile struct from the extraEP dictionary containing full radial data. Sets the multi-radial EP-specific fields plus geometry/physics fields; species AS/TAUS are normalized to electrons (ni/ne and Ti/Te).

preprocess_gacode_inputs(gacode_file, tglfep_file)

Build in-memory Options and profile from input.gacode + input.TGLFEP. No dump.gacode, input.MTGLF, or input.EXPRO required.

Returns (Options, profile, expro_state) where expro_state holds EP vectors for α post-processing and F_REAL setup.

preprocess_imas_inputs(dd, rho, OptionsDict; kw...)

Build TJLFEP's Options/profile inputs directly from an IMAS data dictionary (dd) on the requested rho grid, mirroring the input.gacode preprocessing.

Method provided by the TJLFEPIMASExt package extension, which loads only when IMAS, GACODE, and TurbulentTransport are all present (e.g. under FUSE); calling it without the extension raises a MethodError.

profile{T<:Real}

Local plasma-equilibrium profiles on the TJLFEP radial grid — the intermediary between the input readers and the per-radius TGLF/TJLF solves.

Stores the field signs (SIGN_BT, SIGN_IT), grid/species sizes (NR, NS, N_ION), the per-species matrices (ZS, MASS, AS, TAUS, RLNS, RLTS, VPAR, VPAR_SHEAR), the per-radius geometry/Miller quantities (RMIN, RMAJ, Q, SHEAR, Q_PRIME, P_PRIME, KAPPA/DELTA/ZETA and their shears, ...), and derived quantities (ZEFF, BETAE, RHO_STAR, OMEGA_TAE, gammaE, B_UNIT, ...). Populated from an input.gacode/EXPRO or input.profile source by the readers (e.g. read_input_profile, preprocess_gacode_inputs).

profile_from_gacode(gacode_file; is_ep=2, tglfep_nion=2, q_scale=1.0, rotation_flag=0)

Build a profile struct from input.gacode / dump.gacode, mirroring Fortran expro_read + expro_compute_derived + TGLFEP_read_EXPRO.f90 (species, geometry, gradients). Intended for INPUT_PROFILE_METHOD == 2 parity with Fortran TGLF-EP.

ql_flux_scan_at_marginal(inputsEP, inputsPR; factors, use_gpu) -> (ep_flux, dep, factors)

Run TJLFEP_ky at fixed marginal KYMARK/WIDTH_IN for each scalefactor n_sf. MODE_IN=2 (EP drive only, thermal a/Ln → 1e-6) matches Fortran “thermal gradients off”.

readEXPRO function is a temporary function just used to define any EXPRO constants that are needed. is_EP here is the EXPRO species index (1-based, including electron as species 1 in file layout). For GACODE energetic-ion index use read_expro_for_alpha instead.

readMTGLF: Directly converts a known MTGLF file into a profile struct. It used to do this directly to an inputTJLF struct, but I am going to follow more directly what TGLFEP does first and then come back and make alterations later if I so decide.

Inputs: filepath of input.MTGLF

Outputs: InputTJLF struct to be run through TJLF (as is done in TGLFEP)

readTGLFEP extracts the values needed from the input.TGLFEP file

Inputs: filename

Outputs: Options struct

read_crit_grad(filepath; code="auto") -> (header, values)

Read a TGLF-EP critical-gradient profile file (alpha_dndr_crit.input / alpha_dpdr_crit.input).

Two on-disk layouts are supported and auto-detected (code="auto", the default):

• Fortran TGLFEP/Alpha format — one value per line (F12.4) after the header (this is what TJLFEP now writes; see write_crit_grad).
• Legacy Julia array-literal format — line 2 is [v1, v2, ...].

Pass code="julia" or code="fortran" to force a specific layout.

read_expro_for_alpha(expro_file, prof, is_EP_gacode; gacode_file=nothing)

EXPRO vectors for α post-processing matching Fortran TGLFEP_read_EXPRO:

• GACODE is_EP → fast-ion column in dump.gacode (ni, ti in 10^19/m^3 and keV).
• dlnnidr, dlntidr from expro_bound_deriv(-log(n), rmin) (not profile RLNS).
• max(dlnnidr, 1.0) on the EP species only.

read_gacode_ion_field(path, field, ion_index, nr) -> Vector{Float64}

Read one radial column from a GACODE dump.gacode block (ni, ti, …). ion_index is 1..nion (GACODE ion index, not EXPRO species).

read_gacode_scalar_field(path, field, nr; default=nothing) -> Vector{Float64}

Read a single-column GACODE dump.gacode block (ne, te, rmin, …).

If the block is absent and default is given (not nothing), return fill(default, nr) instead of erroring. Used for optional fields like w0 (toroidal rotation), which some input.gacode files omit; with rotation_flag = 0 a missing w0 is equivalent to zero rotation (it only feeds gammaE/gammap, unused unless ROTATIONAL_SUPPRESSION_FLAG == 1).

read_input_profile(filename::String) -> profile

Read Fortran/GACODE input.profile or dump.profile (fixed columns, NR radial points).

readline_values(file, n_count) -> (v1, v2)

Read one line from an open input.TJLFEP.generate file and parse it as either a scalar or a (min, max) scan range. A leading ( denotes a (v1, v2) range (returned as-is); otherwise the single value is returned duplicated as (val, val). Used by TJLFEP_generate_input.

(n_count is the running count of range-valued parameters; it is incremented locally when a range is parsed.)

readprofile(filename::String) -> profile

Alias for read_input_profile: read an input.profile-style file into a profile struct.

remap_extraEP_for_fortran_save!(args...; kw...)

Remap the extra energetic-particle species bookkeeping into the layout the Fortran-style file save expects (used when persisting IMAS-derived inputs).

Method provided by the TJLFEPIMASExt package extension (loaded under FUSE/IMAS); a MethodError is raised if the extension is not loaded.

runTHD(tglfepfilepath, mtglffilepath, exprofilepath; kwargs...)

Run a TJLFEP critical-gradient scan from the file-based inputs: an input.TGLFEP (scan control), an input.MTGLF (profiles), and an input.EXPRO (experimental profiles for the α/EP quantities).

Reads the three files into Options/profile, applies the EXPRO-derived setup, and runs the scan.

Keywords

• printout::Bool=false - also write an out.TGLFEP summary file
• use_gpu::Bool=false - offload eigensolves to CUDA (auto-detected via pick_device)
• parallel::Symbol=:auto - radial parallel layout
• inner::Symbol=:threads - per-radius backend (:threads or :mps_team)
• team - GPU ids for the MPS team layout
• ql_flux_scan::Bool - also extract QL fluxes at the marginal point
• solver::Symbol=:grid - critical-factor solver (:grid, :ad, :truth, ...)
• refine_rounds::Int=1 - AD (ky,width) refinement passes
• nmodes::Int - number of eigenmodes kept per solve (default from TJLFEP_NMODES, else 4). Applies to all solvers; raise it to catch a higher (e.g. 5th) AE mode.
• k_max::Int - number of :grid scale-factor bracketing rounds (default from TJLFEP_K_MAX, else 4). :grid solver only (inert for the AD solvers).

Returns (width, kymark_out, SFmin, dpdr_crit_out, dndr_crit_out).

An IMAS data-dictionary method runTHD(dd, rho, OptionsDict; ...) is also available when the TJLFEPIMASExt extension is loaded (e.g. under FUSE).

runTHD_dd_radius(args...; kw...)

Single-program-multiple-data (SPMD) per-radius entry point used by the MPS-team multi-GPU layout: runs one radial point of an IMAS-dd scan on a worker.

Method provided by the TJLFEPIMASExt package extension (loaded under FUSE/IMAS); a MethodError is raised if the extension is not loaded.

runTHD_from_gacode(gacode_file, tglfep_file; kwargs...)

Run TJLFEP from input.gacode + input.TGLFEP only. No dump.gacode, input.MTGLF, or input.EXPRO required.

Example

using TJLFEP
width, kymark, SFmin, dpdr, dndr = runTHD_from_gacode(
    "case/input.gacode",
    "case/input.TGLFEP";
    printout=true,
)

run_gacode_scan_task(gacode_file, tglfep_file, scan_index; out_dir, use_gpu, printout)

Run one radius of a SCAN_N job (for Slurm --array=0-(SCAN_N-1)). Writes task_<scan_index>.jls under out_dir with sfmin, width, kymark, and optional out.scalefactor_r### when printout=true.

Use finalize_gacode_scan after all array tasks finish to build α profiles.

save_EXPRO(extraEP::Dict, filename::AbstractString)

Write the extraEP dictionary to a file in input.EXPRO format.

Two-index (species + radial) fields — ni, Ti, dlnnidr, dlntidr: EXPRO_<field>_ <is>_ <ir>= <value>

Single-index (radial) fields — cs, rmin, gammaE, gammap, omegaGAM: EXPRO_<field>_ <ir>= <value>

save_MTGLF(prof::profile, ir_exp::Vector{Int}, filename::AbstractString)

Write a profile struct to a file in input.MTGLF format.

• Scalars (Int, Float): FIELD_NAME= <value>
• Per-species vectors (length NS): FIELD_ <is>= <value>
• Radial vectors (length NR): FIELD_ <ir>= <value>
• Matrices (NR×NS): FIELD_ <ir>_<is>= <value>
• IR_EXP entries appended at end, one per scan radius.

save_TGLFEP(opts::Options, filename::AbstractString)

Write an Options struct to a file in input.TGLFEP format. Format per line: <value> <FIELD_NAME> Scalars (Int, Float, Bool) are written; Vector fields and missing values are skipped, as TGLFEP input files do not encode WIDTH/FACTOR scan vectors.

save_all(opts::Options, prof::profile, extraEP::Dict, dir::AbstractString)

Write all three input files (input.TGLFEP, input.MTGLF, input.EXPRO) to dir.

save_imas_preprocessed_inputs(args...; kw...)

Write the IMAS-derived TJLFEP inputs to disk (the file-based input.TGLFEP / input.MTGLF / input.EXPRO set) so an IMAS run can be reproduced through the file path.

Method provided by the TJLFEPIMASExt package extension (loaded under FUSE/IMAS); a MethodError is raised if the extension is not loaded.

Write input.MTGLF, input.EXPRO from Fortran dump.profile + input.TGLFEP.

setup_gacode_file_inputs(gacode_file, out_dir; tglfep_file, is_ep, ...)

Like setup_fortran_file_inputs, but build the profile from input.gacode instead of a pre-generated dump.profile.

Slurm task id (0-based): SLURM_ARRAY_TASK_ID for --array jobs, else SLURM_PROCID for multi-task jobs (srun -n SCAN_N), else 0.

tjlfepcompleteoutput

inputs: profile_in, inputsEP, inputsPR

outputs: profilein, profileout, irmin, irmax, lacceptprofile

tjlfepcompleteoutput performs the interpolation of a specified profile given a set of radial locations tested. profile_in takes a vector of total ir values and outputs a 201 long vector with this information.

write_crit_grad(io::IO, header, values)
write_crit_grad(path::AbstractString, header, values)

Write a TGLF-EP critical-gradient profile (alpha_dndr_crit.input / alpha_dpdr_crit.input) in the Fortran TGLFEP_driver layout so the Fortran Alpha solver can read it: a one-line header, then one value per line formatted as F12.4.

This matches TGLFEP_driver.f90 (write(22,'(F12.4)') <profile>) and is what Alpha_read_input.f90 expects (it reads the first 5 header chars — Densi / Press — followed by read(6,'(F12.4,A6)') per value). The previous TJLFEP behaviour (println(io, values), a single-line Julia array literal) is not readable by the Fortran Alpha; use this helper instead.