Equilibrium - FEATURE - Accept supplied Coulomb logarithm and collision rates in kinetic profiles

docs/src/conventions.md

@@ -343,9 +343,13 @@ The `KineticForces` module (the Julia port of PENTRC) conventions below are used
 - ``\omega_E`` is the input E × B toroidal-rotation profile, evaluated as
   `welec = omegaE_spline(psi)` (`Torque.jl`). Positive ``\omega_E`` means rotation in the direction
   of the toroidal coordinate ``\zeta``.
-- The mapping to current direction follows from the ``\zeta``/``\phi`` convention fixed by helicity
-  ([Helicity and Handedness](@ref)); it is not re-derived inside the kinetic module, which only
-  reads the supplied profile.
+- In the laboratory frame, where ``\phi`` increases counter-clockwise viewed from above, the toroidal
+  coordinate maps as ``\phi = -h\,(2\pi\zeta + \nu)`` with ``h`` the helicity
+  ([Helicity and Handedness](@ref)), so positive ``\omega_E`` has counter-clockwise sign ``-h``.
+  Relative to the plasma current its sign is ``-h\,s_I = -s_B``, where ``s_I`` and ``s_B`` are the
+  signs of ``I_p`` and ``B_t`` (`ip_sign`, `bt_sign`): positive ``\omega_E`` is co-current exactly when
+  ``B_t`` is clockwise (``s_B = -1``), whatever the helicity. For ``s_I = s_B = +1`` it is clockwise and
+  counter-current.
 
 ### Diamagnetic Frequencies
 

docs/src/kinetic_forces.md

@@ -14,8 +14,9 @@ extension:
     `T_e`, `omega_E` (and `n_i`, defaulting to `n_e` if omitted); optional
     `omega_tor`, `chi_e` and `chi_i` (electron and ion perpendicular heat diffusivities
     ``\chi_{\perp,e}``, ``\chi_{\perp,i}``), `chi_phi` (toroidal momentum diffusivity
-    ``\chi_\phi``), and `zeff`
-    (effective charge ``Z_\mathrm{eff}``, dimensionless). Each dataset
+    ``\chi_\phi``), `zeff`
+    (effective charge ``Z_\mathrm{eff}``, dimensionless), and the supplied collision
+    rates `lnLambda` (dimensionless), `nu_i` and `nu_e` (1/s). Each dataset
     carries a `units` attribute; the root carries `schema_version` and
     `provenance`. Densities are m⁻³, temperatures eV, frequencies rad/s,
     diffusivities m²/s. Additional datasets named `n_*` (e.g. `n_D`, `n_T`)
@@ -24,7 +25,21 @@ extension:
     `write_kinetic_h5` (which round-trips the per-species densities).
   - **ASCII (`.gpeckf`/`.kin`/`.dat`)** — legacy six-column whitespace table
     `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, retained for backward
-    compatibility. Header rows are skipped.
+    compatibility, or, with `kinetic_format = "rates-v1"` in `[Equilibrium]`, nine columns that append
+    `lnLambda nu_i[1/s] nu_e[1/s]` (the Fortran PENTRC rates-v1 table). The column count must match
+    `kinetic_format`. Header rows are skipped.
+
+Supplied collision rates replace the internal Krook collision model of the NTV calculation
+(Logan & Park 2013 Eq. 6), which otherwise computes them from ``n`` and ``T``:
+
+  - Each of `lnLambda`, `nu_i` and `nu_e` is optional and must be positive.
+  - They are interpolated in log space. Outside the table's ψ range the end interval is extended in log
+    space, so supply them over the whole ψ range.
+  - A supplied `nu_i` is the first main ion's rate; the other ion species and the impurity keep the
+    model's ``(Z_s/Z_1)^2\sqrt{m_1/m_s}`` ratio to it.
+  - They feed only the NTV calculation; the tearing layers' resistivity keeps its `lnLambda_form` closure.
+  - The density and temperature scaling factors of `[KineticForces]` do not rescale supplied rates;
+    `collisionality_factor` still multiplies them.
 
 The NTV calculation consumes `n_i, n_e, T_i, T_e, omega_E`; `chi_e`/`chi_i`/`chi_phi`/`zeff`,
 when present, are carried for the resistive-layer (SLAYER) analysis and ignored

docs/src/perturbed_equilibrium.md

@@ -70,10 +70,48 @@ convention the forcing loaders and coil integration produce and the b̃ basis th
 `coupling_overlap` takes care of the conjugation. Singular vectors carry an arbitrary global
 phase, so compare `abs` of overlaps across runs, not the complex value.
 
+In a multi-n run the matrix is block-diagonal in n, so every singular mode has a single n. The ranking spans all n,
+so the dominant mode belongs to the most strongly coupled n, and `DominantMode` stores that mode. For one n, decompose
+its block: `dominant_coupling(rc.C[rows, cols], rc.rational_psi[rows])` with `rows = rc.rational_n .== n` and
+`cols = rc.n_modes .== n`, and project `rootarea_field(rc, modes)[cols]`.
+
 As a summary the run also stores the core-window decomposition (`ψ_N ≤ CORE_PSI_HIGH = 0.9`) under
 `PerturbedEquilibrium/SingularCoupling/DominantMode/`, with `forcing_overlap` holding the run's
 own forcing coefficients `Vᴴ·b̃_x`.
 
+### Other coupling matrices and local windows
+
+The run writes up to five coupling matrices under `PerturbedEquilibrium/SingularCoupling/`. Each has one row per
+rational surface (labelled by `rational_psi`, `rational_q`, `rational_m`, `rational_n`) and acts on the same
+applied b̃ columns (`C_penetrated_area_weighted_field` only when the run computes penetrated fields):
+
+| Dataset | Maps b̃ to | Units |
+|:--- |:--- |:--- |
+| `C_resonant_area_weighted_field` | resonant area-weighted field b̄ʳ = Φʳ/Aʳ | dimensionless |
+| `C_resonant_current` | pitch-resonant current | A/T |
+| `C_island_width_sq` | squared island half-width | 1/T |
+| `C_penetrated_area_weighted_field` | penetrated area-weighted field | dimensionless |
+| `C_Delta_prime` | forcing-driven Δ′ | 1/T |
+
+`ResonantCoupling` wraps the first. Any of the others decomposes the same way through the bare-matrix form
+of `dominant_coupling`; take the b̃ conform from a `ResonantCoupling` of the same file:
+
+```julia
+using HDF5
+C, ψr = h5open("gpec.h5") do f
+    g = f["PerturbedEquilibrium/SingularCoupling"]
+    read(g["C_resonant_current"]), read(g["rational_psi"])
+end
+dom = dominant_coupling(C, ψr; psi_low=0.3, psi_high=0.6)  # any window of surfaces
+c = coupling_overlap(dom, rootarea_field(rc, coil_modes))   # rc = ResonantCoupling("gpec.h5")
+dom.singular_values .* c                                    # the quantity each mode drives, in C's units times T
+```
+
+The bare form carries no (m, n) labels, so `check_mode_basis` cannot guard the projection: build `rc` from
+the same file as `C`.
+
+Any window, on any of the five matrices, is one SVD of the stored matrix.
+
 ```@docs
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.ResonantCoupling
 GeneralizedPerturbedEquilibrium.PerturbedEquilibrium.DominantCoupling

examples/TJ_kinetic_NTV_example/kinetic_rates.dat

@@ -0,0 +1,44 @@
+# kinetic.dat with supplied collision rates (kinetic_format = "rates-v1"): lnLambda is the internal
+# collision model's, and nu_i, nu_e are 3x the model's, so a run that ignores them moves its NTV torque.
+# psi_n n_i n_e T_i[eV] T_e[eV] omega_E lnLambda nu_i[1/s] nu_e[1/s]
+0.000000 1.223086e+19 1.223086e+19 200.0000 200.0000 20000.0000 15.936447 1.320814e+04 8.004180e+05
+0.025000 1.222817e+19 1.222817e+19 195.5000 195.5000 19500.0000 15.902471 1.363412e+04 8.262323e+05
+0.050000 1.203522e+19 1.203522e+19 191.0000 191.0000 19000.0000 15.875444 1.387293e+04 8.407042e+05
+0.075000 1.183797e+19 1.183797e+19 186.5000 186.5000 18500.0000 15.847943 1.411790e+04 8.555497e+05
+0.100000 1.163625e+19 1.163625e+19 182.0000 182.0000 18000.0000 15.819900 1.436971e+04 8.708092e+05
+0.125000 1.142988e+19 1.142988e+19 177.5000 177.5000 17500.0000 15.791293 1.462851e+04 8.864927e+05
+0.150000 1.121868e+19 1.121868e+19 173.0000 173.0000 17000.0000 15.762099 1.489446e+04 9.026097e+05
+0.175000 1.100245e+19 1.100245e+19 168.5000 168.5000 16500.0000 15.732297 1.516771e+04 9.191682e+05
+0.200000 1.078097e+19 1.078097e+19 164.0000 164.0000 16000.0000 15.701861 1.544832e+04 9.361739e+05
+0.225000 1.055403e+19 1.055403e+19 159.5000 159.5000 15500.0000 15.670764 1.573641e+04 9.536320e+05
+0.250000 1.032142e+19 1.032142e+19 155.0000 155.0000 15000.0000 15.638979 1.603203e+04 9.715466e+05
+0.275000 1.008285e+19 1.008285e+19 150.5000 150.5000 14500.0000 15.606479 1.633510e+04 9.899126e+05
+0.300000 9.838102e+18 9.838102e+18 146.0000 146.0000 14000.0000 15.573231 1.664558e+04 1.008728e+06
+0.325000 9.586866e+18 9.586866e+18 141.5000 141.5000 13500.0000 15.539205 1.696325e+04 1.027979e+06
+0.350000 9.328821e+18 9.328821e+18 137.0000 137.0000 13000.0000 15.504369 1.728774e+04 1.047643e+06
+0.375000 9.063714e+18 9.063714e+18 132.5000 132.5000 12500.0000 15.468687 1.761871e+04 1.067700e+06
+0.400000 8.791260e+18 8.791260e+18 128.0000 128.0000 12000.0000 15.432119 1.795560e+04 1.088116e+06
+0.425000 8.511117e+18 8.511117e+18 123.5000 123.5000 11500.0000 15.394627 1.829757e+04 1.108839e+06
+0.450000 8.222884e+18 8.222884e+18 119.0000 119.0000 11000.0000 15.356177 1.864340e+04 1.129796e+06
+0.475000 7.926125e+18 7.926125e+18 114.5000 114.5000 10500.0000 15.316732 1.899140e+04 1.150886e+06
+0.500000 7.620476e+18 7.620476e+18 110.0000 110.0000 10000.0000 15.276253 1.933963e+04 1.171988e+06
+0.525000 7.305553e+18 7.305553e+18 105.5000 105.5000 9500.0000 15.234701 1.968551e+04 1.192949e+06
+0.550000 6.980924e+18 6.980924e+18 101.0000 101.0000 9000.0000 15.192042 2.002559e+04 1.213557e+06
+0.575000 6.646101e+18 6.646101e+18 96.5000 96.5000 8500.0000 15.148252 2.035526e+04 1.233536e+06
+0.600000 6.300621e+18 6.300621e+18 92.0000 92.0000 8000.0000 15.103311 2.066864e+04 1.252527e+06
+0.625000 5.944099e+18 5.944099e+18 87.5000 87.5000 7500.0000 15.057211 2.095834e+04 1.270082e+06
+0.650000 5.576149e+18 5.576149e+18 83.0000 83.0000 7000.0000 15.009964 2.121461e+04 1.285613e+06
+0.675000 5.196380e+18 5.196380e+18 78.5000 78.5000 6500.0000 14.961619 2.142465e+04 1.298341e+06
+0.700000 4.804420e+18 4.804420e+18 74.0000 74.0000 6000.0000 14.912282 2.157127e+04 1.307226e+06
+0.725000 4.399955e+18 4.399955e+18 69.5000 69.5000 5500.0000 14.862145 2.163009e+04 1.310791e+06
+0.750000 3.982956e+18 3.982956e+18 65.0000 65.0000 5000.0000 14.811521 2.157612e+04 1.307520e+06
+0.775000 3.553750e+18 3.553750e+18 60.5000 60.5000 4500.0000 14.760915 2.136511e+04 1.294733e+06
+0.800000 3.113088e+18 3.113088e+18 56.0000 56.0000 4000.0000 14.711172 2.094570e+04 1.269317e+06
+0.825000 2.662383e+18 2.662383e+18 51.5000 51.5000 3500.0000 14.663713 2.024605e+04 1.226918e+06
+0.850000 2.204114e+18 2.204114e+18 47.0000 47.0000 3000.0000 14.621011 1.916915e+04 1.161657e+06
+0.875000 1.742718e+18 1.742718e+18 42.5000 42.5000 2500.0000 14.587459 1.758553e+04 1.065689e+06
+0.900000 1.285949e+18 1.285949e+18 38.0000 38.0000 2000.0000 14.571583 1.533181e+04 9.291128e+05
+0.925000 8.470999e+17 8.470999e+17 33.5000 33.5000 1500.0000 14.591339 1.221607e+04 7.402979e+05
+0.950000 4.497124e+17 4.497124e+17 29.0000 29.0000 1000.0000 14.691469 8.108515e+03 4.913788e+05
+0.975000 1.378556e+17 1.378556e+17 24.5000 24.5000 500.0000 15.024582 3.302452e+03 2.001298e+05
+1.000000 2.762506e+16 2.762506e+16 20.0000 20.0000 0.0000 15.529067 9.192673e+02 5.570792e+04

regression-harness/cases/tj_kinetic_ntv_rates.toml

@@ -0,0 +1,40 @@
+# Regression case: the TJ-analytic n=1 NTV torque with supplied collision rates.
+# Reuses the TJ kinetic NTV deck via [overrides], reading kinetic_rates.dat as a rates-v1 table: its lnΛ is the
+# internal model's and its ν_i, ν_e are 3× the model's, so a reader or spline path that drops the supplied rates
+# moves the torque away from this case's reference. The only harness coverage of supplied collision rates.
+# Each [quantities.*] block names an HDF5 path in the run output, how to extract it, and the
+# noise floor below which a difference is treated as zero.
+[case]
+name = "tj_kinetic_ntv_rates"
+description = "TJ-analytic equilibrium, n=1, no wall, NTV torque with supplied rates-v1 collision rates (reuses the TJ NTV deck via overrides)"
+example_dir = "examples/TJ_kinetic_NTV_example"
+
+[overrides]
+"Equilibrium.kinetic_file" = "kinetic_rates.dat"   # the same profiles with supplied lnLambda, nu_i, nu_e
+"Equilibrium.kinetic_format" = "rates-v1"          # read columns 7-9 as supplied collision rates
+
+# NTV torque with the supplied rates. [Re, Im] pair; deterministic for fixed code and inputs, so the
+# noise floor sits near FP noise.
+[quantities.ntv_torque]
+h5path = "KineticForces/fgar/total_torque"
+type = "real_vector"
+extract = "all_real"
+label = "NTV torque fgar with supplied rates [Re, Im]"
+noise_threshold = 1e-8
+order = 10
+
+[quantities.ntv_psi_nsteps]
+h5path = "KineticForces/fgar/psi_nsteps"
+type = "int_scalar"
+extract = "value"
+label = "NTV ψ quadrature evaluations"
+noise_threshold = 0
+order = 11
+
+# Runtime (special: not from H5)
+[quantities.runtime]
+h5path = ""
+type = "runtime"
+extract = "value"
+label = "Runtime (s)"
+noise_threshold = 0.0

src/Equilibrium/Config.jl

@@ -60,6 +60,8 @@ specified in the input.
   - `kinetic_file::String` - Kinetic-profile file read once into `PlasmaEquilibrium.kinetic` (see [`read_kinetic_file`](@ref));
     empty for none. Every kinetic consumer (NTV, tearing, the layer-overlap scan, grid refinement) derives its view from it.
   - `kinetic_group::String` - HDF5 group of `kinetic_file` holding the profiles
+  - `kinetic_format::String` - Column layout of an ASCII `kinetic_file`: `"legacy"` (6 columns) or `"rates-v1"` (9 columns,
+    appending supplied `lnLambda nu_i nu_e`); see [`read_kinetic_file`](@ref)
   - `ion_species`, `impurity_charge`, `impurity_mass`, `zeff_default`, `chi_e_default`, `chi_i_default`,
     `chi_phi_default`, `tau_E`, `resistivity_model`, `lnLambda_form` - The plasma composition, profile defaults
     and resistivity closure applied to the profiles of `kinetic_file` (see [`KineticProfileData`](@ref)); `nothing`
@@ -96,6 +98,7 @@ specified in the input.
 
     kinetic_file::String = ""
     kinetic_group::String = "/"
+    kinetic_format::String = "legacy"
     ion_species::Union{Nothing,Vector{IonSpecies}} = nothing
     impurity_charge::Union{Nothing,Int} = nothing
     impurity_mass::Union{Nothing,Int} = nothing
@@ -114,7 +117,7 @@ specified in the input.
     function EquilibriumConfig(eq_type, eq_filename, r0_norm, b0_norm, jac_type, _, _, _, _,
         jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc,
         grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, flux_surface_rtol,
-        imas_cocos, kinetic_file, kinetic_group, composition...)
+        imas_cocos, kinetic_file, kinetic_group, kinetic_format, composition...)
         if jac_type == "hamada"
             @info "Forcing hamada coordinate jacobian exponents: power_*"
             power_b = 0
@@ -179,15 +182,20 @@ specified in the input.
             @warn "psihigh = $psihigh exceeds 1.0 (separatrix); clamping to 1.0"
         end
         psihigh = min(psihigh, 1.0)
+        kinetic_format in KINETIC_FORMATS || error("kinetic_format = \"$kinetic_format\" must be one of $KINETIC_FORMATS")
         isempty(kinetic_file) && any(!isnothing, composition) &&
             error("the kinetic composition keys of [Equilibrium] describe the profiles of kinetic_file; set kinetic_file")
         return new(eq_type, eq_filename, r0_norm, b0_norm, jac_type, power_bp, power_b, power_r, power_rc,
             jac_custom_power_bp, jac_custom_power_b, jac_custom_power_r, jac_custom_power_rc,
             grid_type, psilow, psihigh, mpsi, psi_accuracy, mtheta, newq0, flux_surface_rtol,
-            imas_cocos, kinetic_file, kinetic_group, composition...)
+            imas_cocos, kinetic_file, kinetic_group, kinetic_format, composition...)
     end
 end
 
+# ASCII kinetic-table layouts, named as Fortran PENTRC `kinetic_format`, and their column counts.
+const KINETIC_FORMATS = ("legacy", "rates-v1")
+const KINETIC_FORMAT_COLUMNS = Dict("legacy" => 6, "rates-v1" => 9)
+
 # The kinetic composition keys of a config, the ones it sets, as `KineticProfileData` keywords.
 const KINETIC_COMPOSITION_KEYS = (:ion_species, :impurity_charge, :impurity_mass, :zeff_default, :chi_e_default, :chi_i_default,
     :chi_phi_default, :tau_E, :resistivity_model, :lnLambda_form)