diff --git a/docs/src/conventions.md b/docs/src/conventions.md index 154afed892..f599fef855 100644 --- a/docs/src/conventions.md +++ b/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 diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 091492b543..5a315ff9a0 100644 --- a/docs/src/kinetic_forces.md +++ b/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 diff --git a/docs/src/perturbed_equilibrium.md b/docs/src/perturbed_equilibrium.md index 016c70d471..e23e3c5626 100644 --- a/docs/src/perturbed_equilibrium.md +++ b/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 diff --git a/examples/TJ_kinetic_NTV_example/kinetic_rates.dat b/examples/TJ_kinetic_NTV_example/kinetic_rates.dat new file mode 100644 index 0000000000..d786af31c7 --- /dev/null +++ b/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 diff --git a/regression-harness/cases/tj_kinetic_ntv_rates.toml b/regression-harness/cases/tj_kinetic_ntv_rates.toml new file mode 100644 index 0000000000..e4885c1df0 --- /dev/null +++ b/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 diff --git a/src/Equilibrium/Config.jl b/src/Equilibrium/Config.jl index 646c1da2ea..9738fae5b1 100644 --- a/src/Equilibrium/Config.jl +++ b/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) diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index dd910df689..b227f876b6 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -28,11 +28,13 @@ _collision_tkev(T) = max(T / _KEV_J, _T_FLOOR_KEV) # NRL-formulary Coulomb logarithm (natural log) at the floored temperature, floored at 1; n_e in m⁻³, T_e in J. _coulomb_log(ne, Te) = max(17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(_collision_tkev(Te)), 1.0) -# Warn once per kinetic source, as Fortran PENTRC, when a cold node engages the temperature or Coulomb-logarithm floor. -function _warn_collision_floors(source, ne, Te, Ti) +# Warn once per kinetic source, as Fortran PENTRC, when a cold node engages the temperature or Coulomb-logarithm floor +# of a model rate (supplied rates in `data` use neither). +function _warn_collision_floors(source, data, ne, Te, Ti) + data.nu_i !== nothing && data.nu_e !== nothing && return nothing any(T -> T / _KEV_J < _T_FLOOR_KEV, Te) | any(T -> T / _KEV_J < _T_FLOOR_KEV, Ti) && @warn "$source: flooring T at 1 meV in the collision frequencies" maxlog = 1 _id = Symbol("T_floor ", source) - any(i -> 17.3 - 0.5 * log(ne[i] / 1.0e20) + 1.5 * log(_collision_tkev(Te[i])) < 1, eachindex(ne)) && + data.lnLambda === nothing && any(i -> 17.3 - 0.5 * log(ne[i] / 1.0e20) + 1.5 * log(_collision_tkev(Te[i])) < 1, eachindex(ne)) && @warn "$source: flooring the Coulomb logarithm at 1 in the collision frequencies" maxlog = 1 _id = Symbol("lnLambda_floor ", source) return nothing end @@ -40,6 +42,9 @@ end # Names kinetic profiles in error messages by their provenance, when they carry one. _describe(data) = isempty(data.provenance) ? "kinetic profiles" : "kinetic profiles ($(data.provenance))" +# Supplied rate on the kinetic grid: pchip in log space, as Fortran rates-v1, so it stays positive between points. +_supplied_rate(psi_in, v, psi_reg) = v === nothing ? nothing : exp.(_pchip_resample(psi_in, log.(v), psi_reg)) + # Hard stop on an unphysical negative density or temperature (bad input or an extrapolation past the table). function _assert_nonneg_profile(source, label, arr, psi_reg) i = findfirst(<(0), arr) @@ -77,29 +82,37 @@ const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") # The optional profile datasets of the HDF5 kinetic schema and their units, one row per `KineticProfileData` # field of the same name; the reader and writer both iterate it, so a new profile is a row here plus its field. const _KINETIC_H5_PROFILES = ("n_i" => "m^-3", "n_e" => "m^-3", "T_i" => "eV", "T_e" => "eV", "omega_E" => "rad/s", - "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s", "zeff" => "dimensionless") + "omega_tor" => "rad/s", "chi_e" => "m^2/s", "chi_i" => "m^2/s", "chi_phi" => "m^2/s", "zeff" => "dimensionless", + "lnLambda" => "dimensionless", "nu_i" => "1/s", "nu_e" => "1/s") """ - read_kinetic_file(path; group="/") -> KineticProfileData + read_kinetic_file(path; group="/", format="legacy") -> KineticProfileData Read a kinetic-profile file, dispatching on extension: HDF5 (`.h5`/`.hdf5`) via -the GPEC kinetic schema, otherwise a 6-column ASCII table (`.gpeckf`/`.kin`/`.dat`: -`psi_n n_i n_e T_i[eV] T_e[eV] omega_E`). Returns the raw columns; downstream +the GPEC kinetic schema, otherwise an ASCII table (`.gpeckf`/`.kin`/`.dat`) whose layout `format` names: +`"legacy"`, 6 columns `psi_n n_i n_e T_i[eV] T_e[eV] omega_E`, or `"rates-v1"`, 9 columns that append +`lnLambda nu_i[1/s] nu_e[1/s]` (the Fortran PENTRC rates-v1 table). Returns the raw columns; downstream consumers validate the fields they require. """ -function read_kinetic_file(path::AbstractString; group::AbstractString="/") +function read_kinetic_file(path::AbstractString; group::AbstractString="/", format::AbstractString="legacy") isfile(path) || error("Kinetic profile file not found: $path") + format in KINETIC_FORMATS || throw(ArgumentError("kinetic format \"$format\" must be one of $KINETIC_FORMATS")) ext = lowercase(splitext(path)[2]) - return ext in _KINETIC_H5_EXTS ? _read_kinetic_h5(path; group=group) : _read_kinetic_ascii(path) + return ext in _KINETIC_H5_EXTS ? _read_kinetic_h5(path; group=group) : _read_kinetic_ascii(path, format) end """ -Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`. +Read an ASCII kinetic table in the layout `format` names into a `KineticProfileData`. """ -function _read_kinetic_ascii(path::AbstractString) - psi, ni, ne, Ti, Te, omegaE = _read_kinetic_table(path) - return KineticProfileData(; psi=psi, n_i=ni, n_e=ne, T_i=Ti, T_e=Te, omega_E=omegaE, - provenance="ASCII 6-column: $(basename(path))") +function _read_kinetic_ascii(path::AbstractString, format::AbstractString) + cols = _read_kinetic_table(path) + ncol = KINETIC_FORMAT_COLUMNS[format] + length(cols) == ncol || + error("kinetic file $path has $(length(cols)) columns, but kinetic_format = \"$format\" needs $ncol " * + "(\"legacy\": 6; \"rates-v1\": 9, appending lnLambda nu_i nu_e)") + rates = format == "rates-v1" ? (; lnLambda=cols[7], nu_i=cols[8], nu_e=cols[9]) : (;) + return KineticProfileData(; psi=cols[1], n_i=cols[2], n_e=cols[3], T_i=cols[4], T_e=cols[5], omega_E=cols[6], rates..., + provenance="ASCII $(length(cols))-column: $(basename(path))") end """ @@ -129,7 +142,7 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") end """ - attach_kinetic_profiles!(equil, path; group="/", kwargs...) -> equil + attach_kinetic_profiles!(equil, path; group="/", format="legacy", kwargs...) -> equil attach_kinetic_profiles!(equil, data::KineticProfileData) -> equil Attach kinetic profiles to `equil` as `equil.kinetic`. Every kinetic consumer (the NTV torque, the @@ -144,13 +157,14 @@ function attach_kinetic_profiles!(equil::PlasmaEquilibrium, data::KineticProfile return equil end -attach_kinetic_profiles!(equil::PlasmaEquilibrium, path::AbstractString; group::AbstractString="/", kwargs...) = - attach_kinetic_profiles!(equil, KineticProfileData(read_kinetic_file(path; group=group); kwargs...)) +attach_kinetic_profiles!(equil::PlasmaEquilibrium, path::AbstractString; group::AbstractString="/", format::AbstractString="legacy", + kwargs...) = + attach_kinetic_profiles!(equil, KineticProfileData(read_kinetic_file(path; group=group, format=format); kwargs...)) # The kinetic profiles a config names with its composition applied, read once at construction; `nothing` when it names none. configured_kinetic_profiles(config::EquilibriumConfig) = isempty(config.kinetic_file) ? nothing : - KineticProfileData(read_kinetic_file(config.kinetic_file; group=config.kinetic_group); kinetic_composition(config)...) + KineticProfileData(read_kinetic_file(config.kinetic_file; group=config.kinetic_group, format=config.kinetic_format); kinetic_composition(config)...) """ write_kinetic_h5(path, data; group="/", schema_version="1.0", provenance=data.provenance) @@ -229,7 +243,8 @@ paths (the KineticForces quadrature and the self-consistent kinetic-matrix build this set. Errors on any negative density. Each descriptor's `profiles` carries that species' density as `ni_spline` and its full-composition -Krook collision frequency (Logan & Park 2013 Eq. 6) as `nui_spline`; the other profiles are shared. +Krook collision frequency (Logan & Park 2013 Eq. 6) as `nui_spline`, or, when `nu_i` is supplied, that rate +scaled by (z_s/z₁)²·√(m₁/m_s); the other profiles are shared. A main ion's density is `fraction · n_i` (`n_i` is the total main-ion density) or the named `n_*` dataset given by `density` (e.g. `"n_D"`). In a mixed fraction/density list the impurity content is set by the file's n_i/n_e deficit. The `zpitch` closure is exact for z=1 main ions and approximate @@ -264,7 +279,7 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect _assert_positive_density(source, "n_i (total)", ni_total, psi_reg) _assert_nonneg_profile(source, "T_i", Ti, psi_reg) _assert_nonneg_profile(source, "T_e", Te, psi_reg) - _warn_collision_floors(source, ne, Te, Ti) + _warn_collision_floors(source, data, ne, Te, Ti) # Fractions are shares of the total n_i: never above 1, and exactly 1 when every species is fraction-specified. fracs = [s.fraction for s in ion_species if !isnan(s.fraction)] @@ -313,15 +328,17 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect n_main = zeros(npts) n_imp = zeros(npts) nue = zeros(npts) + lnL_sup, nui_sup, nue_sup = (_supplied_rate(psi_in, getfield(data, f), psi_reg) for f in (:lnLambda, :nu_i, :nu_e)) for i in 1:npts z, zp, nm, nimp = multi_ion_composition(zs, [ns[si][i] for si in eachindex(ion_species)], ne[i]; zimp=zimp, mimp=mimp) zeff[i] = z zpitch[i] = zp n_main[i] = nm n_imp[i] = nimp - loglam[i] = _coulomb_log(ne[i], Te[i]) + loglam[i] = lnL_sup === nothing ? _coulomb_log(ne[i], Te[i]) : lnL_sup[i] nue[i] = (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * _collision_tkev(Te[i])^1.5) end + nue_sup === nothing || (nue .= nue_sup) # A negative impurity density means the main ions over-neutralize (Σ z_s n_s > n_e) — unphysical. _assert_nonneg_profile(source, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg) # Approaching the zpitch pole at Zeff = zimp signals a dilute main-ion mix taken up by impurity. @@ -330,7 +347,12 @@ function resolve_ntv_species(data::KineticProfileData, ion_species::AbstractVect # Per-species Krook ν (Logan & Park 2013 Eq. 6) with a test-particle z_s²; `zpitch` is a main-ion closure, # applied approximately to the impurity/electron test species. - _nu(zsp, msp) = [(zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * _collision_tkev(Ti[i])^1.5) for i in 1:npts] + # A supplied ν_i is the first main ion's rate; the others keep the model's (z_s/z₁)²·√(m₁/m_s) ratio to it. + z1, m1 = ion_species[1].z, ion_species[1].m + _nu(zsp, msp) = + nui_sup === nothing ? + [(zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * _collision_tkev(Ti[i])^1.5) for i in 1:npts] : + nui_sup .* ((zsp / z1)^2 * sqrt(m1 / msp)) _view(dens, nu) = KineticProfileSplines(psi_reg, dens, ne, Ti, Te, omegaE, loglam, nu, nue, zeff) # Full NTV species set: main ions, the neutrality-closing impurity, and (optionally) electrons. @@ -386,8 +408,9 @@ The NTV view of raw kinetic profiles (e.g. `equil.kinetic`): interpolate onto a 101-point ψ grid, optionally apply profile scaling knobs, derive collisional / Z_eff diagnostics, and return a `KineticProfileSplines` of named splines (pchip for lnΛ and ν, cubic otherwise). -Needs `n_e`, `T_i`, `T_e` and `omega_E`; `n_i` defaults to `n_e`. See -[`read_kinetic_file`](@ref) for the file formats. +Needs `n_e`, `T_i`, `T_e` and `omega_E`; `n_i` defaults to `n_e`. Supplied `lnLambda`, `nu_i` and `nu_e` +replace the model's values and are not rescaled by the scaling factors. See [`read_kinetic_file`](@ref) +for the file formats. # Arguments @@ -482,12 +505,13 @@ function kinetic_profile_splines(data::KineticProfileData; omegaE .*= ExB_rotation_factor end # After scaling, so a small temperature_factor that engages the floor is reported. - _warn_collision_floors(_describe(data), ne, Te, Ti) + _warn_collision_floors(_describe(data), data, ne, Te, Ti) loglam = zeros(Float64, nkin + 1) nui = zeros(Float64, nkin + 1) nue = zeros(Float64, nkin + 1) zeff = zeros(Float64, nkin + 1) + lnL_sup, nui_sup, nue_sup = (_supplied_rate(psi_input, getfield(data, f), psi_reg) for f in (:lnLambda, :nu_i, :nu_e)) for i in 1:(nkin+1) n_i = ni[i] @@ -498,7 +522,7 @@ function kinetic_profile_splines(data::KineticProfileData; z = zimp - (n_i / n_e) * zi * (zimp - zi) zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z) - ll = _coulomb_log(n_e, T_e) + ll = lnL_sup === nothing ? _coulomb_log(n_e, T_e) : lnL_sup[i] loglam[i] = ll # Krook ν (Logan & Park 2013 Eq. 6) with a test-particle z²; the single-species limit of resolve_ntv_species. @@ -506,32 +530,35 @@ function kinetic_profile_splines(data::KineticProfileData; nue[i] = (zpitch / _NU_PREFAC) * n_e * ll / (sqrt(_ME / _MP) * _collision_tkev(T_e)^1.5) zeff[i] = z end + # Supplied rates replace the model's, as Fortran PENTRC `rates-v1`; the profile scaling factors do not rescale them. + nui_sup === nothing || (nui .= nui_sup) + nue_sup === nothing || (nue .= nue_sup) return KineticProfileSplines(psi_reg, ni, ne, Ti, Te, omegaE, loglam, nui, nue, zeff) end """ -Parse a 6-column kinetic profile table, dropping `#` comments and non-numeric header rows. -Returns the six columns. +Parse a 6- or 9-column kinetic profile table, dropping `#` comments and non-numeric header rows. +Returns the columns. """ function _read_kinetic_table(kinetic_file::AbstractString) table = DelimitedFiles.readdlm(kinetic_file; comments=true) + # readdlm pads a short row with empty strings, so a data row (numeric first entry) that is not all numeric is ragged. + ragged = findfirst(row -> row[1] isa Number && !all(x -> isa(x, Number), row), collect(eachrow(table))) + ragged === nothing || + error("kinetic file $kinetic_file: data row $ragged has fewer columns than the others; every data row needs 6 columns, or 9 with lnLambda nu_i nu_e") + # Keep numeric rows and stack them row-wise (reshape(:, 6) is column-major and would scramble columns). numeric_rows = [collect(row) for row in eachrow(table) if all(x -> isa(x, Number), row)] if isempty(numeric_rows) error("No numeric data rows found in kinetic file: $kinetic_file") end - table = reduce(vcat, (reshape(Float64.(r), 1, 6) for r in numeric_rows)) - - psi_input = collect(table[:, 1]) - n_i_input = collect(table[:, 2]) - n_e_input = collect(table[:, 3]) - T_i_input = collect(table[:, 4]) - T_e_input = collect(table[:, 5]) - omega_e_input = collect(table[:, 6]) - - return psi_input, n_i_input, n_e_input, T_i_input, T_e_input, omega_e_input + ncol = length(numeric_rows[1]) + ncol in (6, 9) || + error("kinetic file $kinetic_file: every data row needs 6 columns, or 9 with lnLambda nu_i nu_e") + table = reduce(vcat, (reshape(Float64.(r), 1, ncol) for r in numeric_rows)) + return [collect(table[:, j]) for j in 1:ncol] end """ diff --git a/src/Equilibrium/PlasmaEquilibrium.jl b/src/Equilibrium/PlasmaEquilibrium.jl index dd1737bf66..52d62d72b2 100644 --- a/src/Equilibrium/PlasmaEquilibrium.jl +++ b/src/Equilibrium/PlasmaEquilibrium.jl @@ -146,6 +146,8 @@ layer-overlap scan) reads the composition and closure from here, so they all see | `chi_i` | ion perpendicular heat diffusivity χ⊥,i | | `chi_phi` | toroidal momentum diffusivity χ_φ | | `zeff` | effective charge Z_eff(ψ) supplied directly (optional); the resistive-layer analysis prefers it to its closure | +| `lnLambda` | Coulomb logarithm supplied directly (optional); replaces the internal one in the NTV collision model | +| `nu_i`/`nu_e` | main-ion / electron collision rates in 1/s supplied directly (optional); replace the internal collision model | | `species_densities` | named per-species densities (e.g. `"n_D"`, `"n_T"`) for explicit multi-ion input; `nothing` for ASCII / no extras | | `provenance` | short string recording the source file/format | | `ion_species` | main-ion species, a `Vector{IonSpecies}`; default one deuterium ion carrying all of `n_i` | @@ -174,6 +176,9 @@ struct KineticProfileData chi_i::Union{Nothing,Vector{Float64}} chi_phi::Union{Nothing,Vector{Float64}} zeff::Union{Nothing,Vector{Float64}} + lnLambda::Union{Nothing,Vector{Float64}} + nu_i::Union{Nothing,Vector{Float64}} + nu_e::Union{Nothing,Vector{Float64}} species_densities::Union{Nothing,Dict{String,Vector{Float64}}} provenance::String ion_species::Vector{IonSpecies} @@ -193,7 +198,7 @@ const _LNLAMBDA_FORMS = (:nrl, :sauter, :wesson) function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=nothing, omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_i=nothing, chi_phi=nothing, zeff=nothing, - species_densities=nothing, provenance="", + lnLambda=nothing, nu_i=nothing, nu_e=nothing, species_densities=nothing, provenance="", ion_species=[IonSpecies(; z=1, m=2, fraction=1.0)], impurity_charge=6, impurity_mass=12, zeff_default=1.0, chi_e_default=1.0, chi_i_default=nothing, chi_phi_default=1.0, tau_E=nothing, resistivity_model=:sauter, lnLambda_form=:nrl) @@ -209,8 +214,16 @@ function KineticProfileData(; psi, n_i=nothing, n_e=nothing, T_i=nothing, T_e=no model, form = Symbol(resistivity_model), Symbol(lnLambda_form) model in _RESISTIVITY_MODELS || throw(ArgumentError("KineticProfileData: resistivity_model = $model not in $_RESISTIVITY_MODELS")) form in _LNLAMBDA_FORMS || throw(ArgumentError("KineticProfileData: lnLambda_form = $form not in $_LNLAMBDA_FORMS")) + # An all-zero rate dataset counts as not supplied, like the other optional profiles. + _rate(x) = x === nothing || all(iszero, x) ? nothing : x + lnLambda, nu_i, nu_e = _rate(lnLambda), _rate(nu_i), _rate(nu_e) + for (name, v) in (("lnLambda", lnLambda), ("nu_i", nu_i), ("nu_e", nu_e)) + v === nothing || all(x -> isfinite(x) && x > 0, v) || + throw(ArgumentError("KineticProfileData: supplied $name must be positive and finite everywhere")) + end return KineticProfileData(Float64.(collect(psi)), _v(n_i), _v(n_e), _v(T_i), _v(T_e), - _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), _v(zeff), species_densities, String(provenance), + _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_i), _v(chi_phi), _v(zeff), _v(lnLambda), _v(nu_i), _v(nu_e), + species_densities, String(provenance), species, Int(impurity_charge), Int(impurity_mass), Float64(zeff_default), _s(chi_e_default), _s(chi_i_default), _s(chi_phi_default), _s(tau_E), model, form) end diff --git a/src/ErrorFields/NTVLimits.jl b/src/ErrorFields/NTVLimits.jl index d75abd21a2..c216f702a5 100644 --- a/src/ErrorFields/NTVLimits.jl +++ b/src/ErrorFields/NTVLimits.jl @@ -530,10 +530,11 @@ function efc_couplings( KineticForces.set_perturbation_data!(kf_intr, pe_state, ffs, equil, ffs.metric) return kf_intr end + # One quiet logger for the whole scan, so a `maxlog` warning prints once per scan rather than once per point. + logger = verbose ? Base.CoreLogging.current_logger() : Base.CoreLogging.SimpleLogger(stderr, Base.CoreLogging.Warn) # Torque and cumulative torque profile (resampled onto the kinetic ψ grid) of one perturbation at one rotation. function torque_of(kf_intr, profiles) kf_state = KineticForces.KineticForcesState() - logger = verbose ? Base.CoreLogging.current_logger() : Base.CoreLogging.SimpleLogger(stderr, Base.CoreLogging.Warn) Base.CoreLogging.with_logger(logger) do KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, profiles) end diff --git a/src/KineticForces/EnergyIntegration.jl b/src/KineticForces/EnergyIntegration.jl index b321f546de..54bf50fdcc 100644 --- a/src/KineticForces/EnergyIntegration.jl +++ b/src/KineticForces/EnergyIntegration.jl @@ -308,7 +308,8 @@ function _integrate_energy_resonant(p::EnergyParams, leff::Float64, wb::Float64, u_max = sqrt(X_ENERGY_MAX) if npole == 0 - val, _ = quadgk(integrand_u, 0.0, u_max; atol=atol, rtol=rtol, segbuf=segbuf) + val, err = quadgk(integrand_u, 0.0, u_max; atol=atol, rtol=rtol, segbuf=segbuf) + check_inner_quadrature(val, err, atol, rtol, "energy") return val end @@ -330,7 +331,8 @@ function _integrate_energy_resonant(p::EnergyParams, leff::Float64, wb::Float64, end sort!(pts) unique!(pts) - smooth_val, _ = _quadgk_breaks(integrand_u, pts, atol, rtol, segbuf) + smooth_val, err = _quadgk_breaks(integrand_u, pts, atol, rtol, segbuf) + check_inner_quadrature(smooth_val, err, atol, rtol, "energy") return smooth_val + pole_contribution end @@ -358,7 +360,8 @@ function integrate_energy(wn::Float64, wt::Float64, we::Float64, wd::Float64, # CGL has no resonance denominator and no pole — integrate the physical # x-space integrand directly over the half line (QuadGK maps [0,∞) itself). if f0type == "cgl" - val, _ = quadgk(x -> _energy_integrand_real(x, p), 0.0, Inf; atol=atol, rtol=rtol, segbuf=segbuf) + val, err = quadgk(x -> _energy_integrand_real(x, p), 0.0, Inf; atol=atol, rtol=rtol, segbuf=segbuf) + check_inner_quadrature(val, err, atol, rtol, "CGL energy") return val end diff --git a/src/KineticForces/PitchIntegration.jl b/src/KineticForces/PitchIntegration.jl index de1714071a..ba09d92773 100644 --- a/src/KineticForces/PitchIntegration.jl +++ b/src/KineticForces/PitchIntegration.jl @@ -96,11 +96,12 @@ function integrate_pitch_gar_quadgk( bobmax_clip = clamp(bobmax, lambda_min, lambda_max) buf = zeros(ComplexF64, nqty) kernel! = (out, λ) -> _pitch_gar_kernel_quadgk!(out, λ, params) - I, _ = if lambda_min < bobmax_clip < lambda_max + I, E = if lambda_min < bobmax_clip < lambda_max quadgk!(kernel!, buf, lambda_min, bobmax_clip, lambda_max; atol=pitch_atol, rtol=pitch_rtol) else quadgk!(kernel!, buf, lambda_min, lambda_max; atol=pitch_atol, rtol=pitch_rtol) end + check_inner_quadrature(I, E, pitch_atol, pitch_rtol, "pitch") return copy(I) end @@ -185,11 +186,12 @@ function integrate_pitch_gar_quadgk_wt( bobmax_clip = clamp(bobmax, lambda_min, lambda_max) buf = zeros(ComplexF64, 2 * nqty) kernel! = (out, λ) -> _pitch_gar_kernel_quadgk_wt!(out, λ, params) - I, _ = if lambda_min < bobmax_clip < lambda_max + I, E = if lambda_min < bobmax_clip < lambda_max quadgk!(kernel!, buf, lambda_min, bobmax_clip, lambda_max; atol=pitch_atol, rtol=pitch_rtol) else quadgk!(kernel!, buf, lambda_min, lambda_max; atol=pitch_atol, rtol=pitch_rtol) end + check_inner_quadrature(I, E, pitch_atol, pitch_rtol, "pitch") return copy(I) end diff --git a/src/KineticForces/Utils.jl b/src/KineticForces/Utils.jl index 60c3aa03e8..5f8a7ae2b8 100644 --- a/src/KineticForces/Utils.jl +++ b/src/KineticForces/Utils.jl @@ -5,6 +5,19 @@ Shared numerical helpers for KineticForces: the generic sign-change root scan an kinetic-resonance surface locator built on it. """ +""" + check_inner_quadrature(I, E, atol, rtol, what) + +Warn once per integral kind `what` and Julia session when the error estimate `E` exceeds `max(atol, rtol·‖I‖)`, +i.e. QuadGK stopped short of its tolerance (evaluation cap or endpoint roundoff), or when `I` or `E` is not finite. +""" +function check_inner_quadrature(I, E, atol, rtol, what::String) + nI = norm(I) + (isfinite(nI) && isfinite(E) && E <= max(atol, rtol * nI)) || + @warn "KineticForces $what integral did not converge: error estimate $E vs |I| = $nI" maxlog = 1 _id = Symbol(what) + return nothing +end + """ find_sign_change_roots(f, grid) → Vector{Float64} diff --git a/test/runtests_kinetic.jl b/test/runtests_kinetic.jl index 2be7e099d7..1053263329 100644 --- a/test/runtests_kinetic.jl +++ b/test/runtests_kinetic.jl @@ -755,6 +755,13 @@ using TOML @test_throws ErrorException KF.check_psi_quadrature_convergence(total, NaN, ctrl, "fgar") end + @testset "check_inner_quadrature" begin + # An error estimate above tolerance (QuadGK stopped short) or a NaN must not pass silently + @test_logs KF.check_inner_quadrature(1.0 + 0im, 1e-6, 0.0, 1e-5, "test converged") + @test_logs (:warn, r"did not converge") KF.check_inner_quadrature(1.0 + 0im, 1e-3, 0.0, 1e-5, "test capped") + @test_logs (:warn, r"did not converge") KF.check_inner_quadrature([1.0 + 0im, NaN], 0.0, 0.0, 1e-5, "test nan") + end + @testset "METHOD_REGISTRY" begin # The registry is the single source of truth for NTV methods. Every entry's # flag must be a real KineticForcesControl field (the TOML kwargs splat relies diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl index 7a1e62c8c6..be3c0e9f5d 100644 --- a/test/runtests_multiion.jl +++ b/test/runtests_multiion.jl @@ -214,6 +214,60 @@ using HDF5 @test minimum(ped_T) <= minimum(resampled) && maximum(resampled) <= maximum(ped_T) end + @testset "supplied collision rates replace the internal model" begin + psi = collect(0.0:0.05:1.0) + prof(a) = @. a * (1 - 0.6psi) + data = EQ.KineticProfileData(; psi=psi, n_i=prof(1.0e20), n_e=prof(1.1e20), T_i=prof(2000.0), T_e=prof(2500.0), + omega_E=prof(1.0e4)) + model = EQ.kinetic_profile_splines(data) + xs = model.xs + # The model's own rates, supplied at the table points, come back to interpolation accuracy. + own = EQ.KineticProfileData(data; lnLambda=model.loglam_spline.(psi), nu_i=model.nui_spline.(psi), nu_e=model.nue_spline.(psi)) + sup = EQ.kinetic_profile_splines(own) + for f in (:loglam_spline, :nui_spline, :nue_spline) + @test getfield(sup, f).(xs) ≈ getfield(model, f).(xs) rtol = 1e-3 + end + # Interpolation is in log space: an exponential in ψ is reproduced exactly, which a linear-space fit would not do. + expo = EQ.kinetic_profile_splines(EQ.KineticProfileData(own; nu_i=@.(1e3 * exp(4psi)), nu_e=@.(1e5 * exp(-3psi)))) + @test expo.nui_spline.(xs) ≈ @.(1e3 * exp(4xs)) rtol = 1e-12 + @test expo.nue_spline.(xs) ≈ @.(1e5 * exp(-3xs)) rtol = 1e-12 + # A supplied lnΛ alone feeds the model's ν, which is linear in it. + lnl = EQ.kinetic_profile_splines(EQ.KineticProfileData(data; lnLambda=fill(15.0, length(psi)))) + @test all(≈(15.0), lnl.loglam_spline.(xs)) + @test lnl.nui_spline.(xs) ≈ model.nui_spline.(xs) .* 15.0 ./ model.loglam_spline.(xs) rtol = 1e-12 + # Multi-ion: the first main ion takes the supplied ν_i; the others keep the model's (z_s/z₁)²·√(m₁/m_s) ratio. + sp = EQ.resolve_ntv_species(own, [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)]; zimp=6, mimp=12) + nu(z, m) = only(s for s in sp if s.z == z && s.m == m).profiles.nui_spline.(xs) + @test nu(1, 2) ≈ sup.nui_spline.(xs) rtol = 1e-12 + @test nu(1, 3) ≈ sqrt(2 / 3) .* nu(1, 2) rtol = 1e-12 + @test nu(6, 12) ≈ 36 * sqrt(2 / 12) .* nu(1, 2) rtol = 1e-12 + # Rates must be positive; an all-zero dataset counts as absent. + @test_throws ArgumentError EQ.KineticProfileData(data; nu_i=-own.nu_i) + @test EQ.KineticProfileData(data; nu_e=zeros(length(psi))).nu_e === nothing + # Both file formats carry the rates: HDF5 round-trips them and a 9-column ASCII table reads them. + h5 = tempname() * ".h5" + EQ.write_kinetic_h5(h5, own) + back = EQ.read_kinetic_file(h5) + @test back.lnLambda ≈ own.lnLambda && back.nu_i ≈ own.nu_i && back.nu_e ≈ own.nu_e + rm(h5) + ascii = tempname() * ".kin" + open(ascii, "w") do io + println(io, "# psi n_i n_e T_i T_e omega_E lnLambda nu_i nu_e") + for k in eachindex(psi) + println(io, join((psi[k], own.n_i[k], own.n_e[k], own.T_i[k], own.T_e[k], own.omega_E[k], own.lnLambda[k], own.nu_i[k], own.nu_e[k]), " ")) + end + end + table = EQ.read_kinetic_file(ascii; format="rates-v1") + @test table.nu_i ≈ own.nu_i && table.lnLambda ≈ own.lnLambda + # The layout comes from the format, not the column count: a 9-column table read as legacy is an error. + @test_throws r"needs 6" EQ.read_kinetic_file(ascii) + @test_throws ArgumentError EQ.read_kinetic_file(ascii; format="rates-v2") + # A data row with fewer columns than the others is an error, not a silently dropped row. + open(io -> println(io, join((1.0, own.n_i[end], own.n_e[end], own.T_i[end], own.T_e[end], own.omega_E[end]), " ")), ascii, "a") + @test_throws r"fewer columns" EQ.read_kinetic_file(ascii; format="rates-v1") + rm(ascii) + end + @testset "combine_species_states" begin # Combined profile = Σ species dT/dψ interpolated onto the union grid; total = Σ total. mk(psi, dt) = ( diff --git a/test/runtests_solve_api.jl b/test/runtests_solve_api.jl index 975a2b57ed..8dd25a951a 100644 --- a/test/runtests_solve_api.jl +++ b/test/runtests_solve_api.jl @@ -386,6 +386,11 @@ using HDF5 @test prob.eta == GPEC.InnerLayer.layer_eta(kgal.equil, kgal.surfaces) && all(>(0), prob.eta) @test prob.rho ≈ fill(2 * GPEC.Utilities.M_P * 5e19, n) @test prob.rotation ≈ fill(1e3, n) + # A supplied Coulomb logarithm feeds only the NTV collision model; η keeps its own closure. + keq = deepcopy(kgal.equil) + attach_kinetic_profiles!(keq, KineticProfileData(; psi, n_e=flat(5e19), n_i=flat(5e19), T_e=flat(1e3), T_i=flat(1e3), + omega_E=flat(2π * 1e3), lnLambda=flat(30.0))) + @test GPEC.InnerLayer.layer_eta(keq, kgal.surfaces) == prob.eta # Explicit values override each default on its own. given = MatchProblem(kgal; eta=1e-7, rotation=fill(3.0, n)) @test given.eta == fill(1e-7, n) && given.rho == prob.rho && given.rotation == fill(3.0, n) @@ -554,6 +559,23 @@ using HDF5 end end + @testset "kinetic_format" begin + # A rates-v1 kinetic file read through the deck's [Equilibrium] keys and through the API keyword. + mktempdir() do dir + psi = range(0.0, 1.0; length=11) + open(joinpath(dir, "rates.kin"), "w") do io + for x in psi + println(io, join((x, 5e19, 5e19, 1e3, 1e3, 1e3, 16.0 + x, 1e3 * exp(x), 1e5 * exp(-x)), " ")) + end + end + eq_overrides = Dict{String,Any}("kinetic_file" => "rates.kin", "kinetic_format" => "rates-v1") + deck_dir = _stage_deck(dir; eq_overrides=eq_overrides) + from_deck = GPEC.Equilibrium.configured_kinetic_profiles(GPEC.equilibrium_config(TOML.parsefile(joinpath(deck_dir, "gpec.toml")), deck_dir)) + from_api = attach_kinetic_profiles!(deepcopy(equil), joinpath(dir, "rates.kin"); format="rates-v1").kinetic + @test from_deck.nu_i == from_api.nu_i && from_deck.lnLambda == from_api.lnLambda && from_deck.nu_e == from_api.nu_e + end + end + @testset "deck kinetic source names" begin # A deck selects only the KineticForces model; test matrices and the retired key stop the parse. parsed(overrides) =