diff --git a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md index 1a3bc1fc..9d17a0bd 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md +++ b/.claude/agent-memory/fortran-physics-reviewer/MEMORY.md @@ -5,3 +5,4 @@ - [KineticForces (NTV) Audit Checklist](kinetic_ntv_map.md) — pentrc->KineticForces map, Logan 2015 matrices, what to verify - [InnerLayer (Resistive) Audit Checklist](resistive_layer_map.md) — rmatch->InnerLayer map, GGJ Wasow basis / Δ′, what to verify - [Galerkin Δ′ Assembly Map](galerkin_assembly_map.md) — gal.f<->GalerkinAssembly.jl; resonant sign chain verified; PASS +- [reg_spot vs singfac_min audit](reg_spot_regularization.md) — reg_spot field-recon smoothing CORRECT; distinct from singfac_min ODE gate diff --git a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md index 414be7a9..a3b73369 100644 --- a/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md +++ b/.claude/agent-memory/fortran-physics-reviewer/kinetic_ntv_map.md @@ -31,3 +31,59 @@ kinetic energy/matrices from trapped-particle nonambipolar transport. Fortran co - **Normalization**: torque Im(T)=2n·δW_k (Eq. 19) and diamagnetic-frequency sign/factor conventions match Fortran `torque.F90`. - **Mode indexing**: m, n ranges and block packing over n stay consistent with ForceFreeStates. - **Method variants**: FGAR/TGAR/PGAR/RLAR/CLAR/FCGL/TMM/WMM each present, not stubbed to a single fallback. + +## Multi-ion (D-T) NTV — composition/collisionality (2026-07) +Both Fortran PENTRC (`inputs.f90:236-243`, `read_kin`) and Julia (`KineticProfiles.jl:261-278`) +support ONE main ion (zi,mi) + ONE impurity (zimp,mimp) per run. Correct multi-main-ion NTV is +an EXTENSION beyond Fortran, but consistent with Logan-Park 2013's pitch-angle (Lorentz) model. +- Zeff = Σ_s Z_s² n_s / n_e; quasineutrality n_e = Σ_s Z_s n_s (all ions incl. impurity). +- Bug in current D-T split (run D, then T, each with ni=Ni/2): the `ni` column sets BOTH species + density AND Zeff via z = zimp-(n_i/n_e)zi(zimp-zi). With ni=Ni/2 it treats the missing half as + high-Z impurity → Zeff≈3.7 instead of true ≈1–1.4. Corrupts zpitch → nue,nui. +- Correct collisionality: ν_a ∝ Z_a² lnΛ · Σ_b n_b Z_b² /(√m_a T_a^{3/2}) = ∝ Z_a² lnΛ·(n_e·Zeff_true). + The zpitch(Zeff) polynomial is a main-ion+impurity closure (momentum-restoring correction); the + MINIMAL fix is to feed the TRUE (full-composition) Zeff into the existing zpitch/ν formulas. +- Additivity: τ = Σ_s τ_s. Lorentz operator is additive over field species; species couple only + through shared δB, shared ω_E, shared Zeff/ν. Additive at this theory's order. +- What changes single→multi: Zeff, zpitch, nue, nui (→nueff). UNCHANGED for equal-shape D-T: + wdian/wdiat (log-derivative, density factor cancels), wtran/wbhat/wdhat/wgyro (already per-species), + and the resonant-density prefactor (Ni/2 per species is correct). +- ASIDE (separate fidelity bug, same code block): Fortran `inputs.f90:238` uses natural log for lnΛ; + Julia `KineticProfiles.jl:270` uses log10 — diverges away from the n=1e20,T=1keV reference point. + +## Multi-ion NTV — full-composition species set (2026-07 audit, PASS-with-caveats) +Reviewed `resolve_ntv_species` (KineticProfiles.jl ~214-266) + `compute_calculated_kinetic_matrices` +(CalculatedKineticMatrices.jl ~96-151). Verdict: physics is sound. +- ν_s field-density RECONCILIATION (supersedes the earlier "should be n_e·Zeff" note above): + code uses `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ/(√m_s·T_i^1.5)`, i.e. field density = zpitch·n_main + (n_main = Σ MAIN-ion densities, no impurity, unweighted). This is CORRECT and MORE faithful to + single-ion PENTRC than n_e·Zeff: PENTRC's design is zpitch·n_i, NOT Σ_b n_b Z_b². Numerically + zpitch·n_main ≈ n_e·Zeff (e.g. Zeff=1.5,C6 → 1.55 vs 1.5); the gap IS the intended momentum- + restoring correction that zpitch(Zeff) carries. So n_main vs n_e·Zeff is immaterial at D-T (z=1). +- z_s² test-particle factor: CORRECT and correctly placed (deflection freq ∝ test charge²). Single-ion + had no z² only because zi=1. Reduces EXACTLY to single-ion nui for one z=1,fraction=1 ion (verified). +- Impurity as its own test species: field density zpitch·n_main is NOT undercounting — the impurity's + z_imp² contribution is already folded into Zeff inside zpitch. CAVEAT: zpitch is strictly a main-ion + momentum-restoring closure; reusing it for the impurity/electron ν is an approximation beyond the + single-ion theory. Acceptable (impurity δW ∝ n_imp is small); worth a one-line annotation. +- Electron descriptor passes `ns[1]` (first ion's density) as its `ni_spline` — HARMLESS: the kernel + `_setup_surface_state` (Torque.jl:659-664) reads `ne_spline` (full shared n_e) when electron=true, + never ni_spline. Electron n_s = full n_e (correct). Cosmetic smell only. +- Self-consistent δW summation (CalculatedKineticMatrices.jl:112-148): kw_flat/kt_flat are PURELY the + kinetic matrices (Logan 7.30-7.35), every term ∝ species phase-space density n_s·f0_s — NO species- + independent baseline. Fluid F,K,G added ONCE downstream in _compute_fkg_matrices!, outside the + species loop. So `+=` over species is clean additivity; NOTHING in kw is wrongly ×species-count. + → TC-24 n=3 δW +0.066(D) → −0.10(D+T+C+e) sign flip is PLAUSIBLE physics near marginal stability, + NOT a double-count. Dominant driver: the newly-ON electron channel (full n_e, opposite precession). + DECISIVE cheap diagnostic to confirm: run D(½)+T(½) with electron=OFF, impurity absent — should + return ≈ +0.066 (single-ion D). If D+T alone ≈ +0.066, the whole shift is electron+impurity = physical. + +## Single-ion nui vs multi-species z_s² (#339, 2026-07 decision) +Fortran PENTRC `inputs.f90:240-241` single-ion nui = (zpitch/3.5e17)·n_i·lnΛ/(√mi·T_i^1.5) has +NO explicit zi² (implicitly assumes zi=1, main ion hydrogenic). Julia `load_kinetic_profiles` +(~L416) is a faithful exact port. Multi-species `_nu` (~L253) adds explicit test-particle z_s² +(the physically correct pitch-angle/Lorentz form; zpitch is the field-side momentum-restoring +factor of Zeff, independent of test charge — no double count). RECOMMENDATION: option (a) ADD +zi² to single-ion nui. It is numerically identical for the default zi=1 (regression byte-identical), +makes single-ion the true 1-species special case of `_nu`, and only "deviates" from Fortran in the +exotic zi≠1 case where Fortran is physically wrong anyway. Annotate as a documented improvement. diff --git a/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md b/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md new file mode 100644 index 00000000..137301e2 --- /dev/null +++ b/.claude/agent-memory/fortran-physics-reviewer/reg_spot_regularization.md @@ -0,0 +1,24 @@ +--- +name: reg_spot vs singfac_min regularization audit +description: Verdict on GPEC-Julia reg_spot (field-reconstruction smoothing) port and its distinction from singfac_min (ODE crossing gate) +metadata: + type: project +--- + +## Verdict (audited FieldReconstruction.jl, no local Fortran repo available) +`reg_spot` port judged CORRECT-WITH-CAVEATS (caveat = Fortran verified from GPEC `xm*` convention + physics + known default, not from on-disk source). + +## reg_spot — field-reconstruction smoothing (PerturbedEquilibrium) +- Factor form: `reg_factor = singfac²/(singfac²+reg_spot²)`, singfac = m - n·q. Correct GPEC form. → 0 at rational surface, → 1 away. Applied in singfac-space, NOT ψ-space. +- Default 5e-2 = 0.05 (PerturbedEquilibriumStructs.jl:50). Matches GPEC gpec_input namelist default. +- Quantities regularized (matches GPEC `xm*` modified-quantity convention): + - xmp1 = ξ^ψ' (FieldReconstruction.jl:393-395) + - xms (clebsch_alpha) computed FROM regularized xmp1 via xms=-A⁻¹(B·xmp1_reg+C·xsp) (line 398-415, gpeq_sol) + - xmt/xmz = regularized ξ^θ,ξ^ζ (line 584-594, gpeq_contra); xwt/xwz kept unregularized + - b^θ_reg,b^ζ_reg inherit reg via xmp1/xms (compute_modified_field_modes) +- LEFT UNREGULARIZED: ξ^ψ primitive (xsp / clebsch_psi = copy, line 356; xwp uses raw xsp). CORRECT — ξ^ψ (normal disp) is finite at rational surface; only 1/singfac-divergent tangential/derivative quantities are softened. + +## singfac_min — DISTINCT mechanism (ForceFreeStates ODE crossing gate) +- Default 1e-4 (ForceFreeStatesStructs.jl:267), "Matches Fortran STRIDE". Gates singular-surface crossings during EL ODE integration (EulerLagrange.jl:466-519, GeneralizedPerturbedEquilibrium.jl:426). +- SingularCoupling.jl:204 uses `spot = 5e-4` as an offset (5e-4/|n·q1|) to evaluate bwp1 at lpsi/rpsi — labeled "matches Fortran default singfac_min". This is the ψ-offset for one-sided derivative eval, a THIRD distinct use, not reg_spot. +- reg_spot (0.05, singfac-space smoothing of reconstructed fields) and singfac_min (1e-4, ODE crossing gate) are separate; both ported. diff --git a/examples/Solovev_kinetic_multiion_example/forcing.dat b/examples/Solovev_kinetic_multiion_example/forcing.dat new file mode 100644 index 00000000..42eb9470 --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/forcing.dat @@ -0,0 +1,5 @@ +# Forcing data for perturbed equilibrium calculations +# normalization: normal_field_T +# Format: n m amplitude_real amplitude_imag +# Single mode test case: n=1, m=2, amplitude=1e-4 T +1 2 1e-4 0.0 diff --git a/examples/Solovev_kinetic_multiion_example/gpec.toml b/examples/Solovev_kinetic_multiion_example/gpec.toml new file mode 100644 index 00000000..b869e04e --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/gpec.toml @@ -0,0 +1,109 @@ +# Solovev analytical equilibrium — n=1 ideal stability + perturbed equilibrium + NTV torque. +# Extends the Solovev ideal example with a [KineticForces] section, so after the plasma +# response is computed the neoclassical toroidal viscosity torque is integrated over ψ as a +# post-processing diagnostic (adaptive quadrature paneled at the rational surfaces). +# The kinetic profiles (kinetic.dat) are constructed to be consistent with this equilibrium's +# pressure; the equilibrium is generated analytically from the embedded [SOL_INPUT] section. + +[Equilibrium] +eq_type = "sol" # Type of the input 2D equilibrium file +jac_type = "pest" # Coordinate system (hamada, pest, boozer, equal_arc, park, other) +grid_type = "ldp" # Radial grid packing type +psilow = 1e-4 # Lower limit of normalized poloidal flux +psihigh = 0.9995 # Upper limit of normalized poloidal flux +mpsi = 128 # Number of radial grid points (0 = auto-compute from psi_accuracy) +mtheta = 256 # Number of poloidal grid points +newq0 = 0 # Override for on-axis safety factor (0 = use input value) +etol = 1e-7 # Error tolerance for equilibrium solver +force_termination = false # Terminate after equilibrium setup (skip stability calculations) + +[Wall] +# Close conformal wall is required to stabilize this Solovev fixture's n=1 external kink: +# with nowall, et[1] = -6.8 (strongly unstable); with this wall, et[1] = +0.24 (barely stable). +# The plasma is near marginal stability, so the BVP Δ' matrix values are pathological +# (dpm magnitudes ~ 10¹¹, |Im/Re| ≫ 1). This fixture's role is integration-pipeline +# smoke testing + et[1] regression, NOT BVP Δ' regression — DIIID-like is the canonical +# Δ'-matrix fixture (stable et[1] = +1.6, clean BVP Δ'). +shape = "conformal" # Wall shape (nowall, conformal, elliptical, dee, mod_dee, filepath) +a = 0.2415 # Distance from plasma (conformal) or shape parameter +aw = 0.05 # Half-thickness parameter for Dee-shaped walls +bw = 1.5 # Elongation parameter for wall shapes +cw = 0 # Offset of wall center from major radius +dw = 0.5 # Triangularity parameter for wall shapes +tw = 0.05 # Sharpness of wall corners (try 0.05 as initial value) +equal_arc_wall = true # Equal arc length distribution of nodes on wall + +[ForcingTerms] +forcing_data_file = "forcing.dat" # Manual mode table (n, m, amplitude_real, amplitude_imag) +forcing_data_format = "ascii" # Format: "ascii", "hdf5", or "coil" (Biot-Savart from 3D wires) + +[PerturbedEquilibrium] +fixed_boundary = false # Use fixed boundary conditions +output_eigenmodes = true # Output eigenmode fields as b-fields +compute_response = true # Compute plasma response to forcing +compute_singular_coupling = true # Compute singular layer coupling metrics +verbose = true # Enable verbose logging +write_outputs_to_HDF5 = true # Write perturbed equilibrium outputs to HDF5 +reg_spot = 0.05 # Regularization width for singular surfaces (0 = disabled) + +[ForceFreeStates] +local_stability_flag = true # Perform local stability analysis (Mercier and ballooning) across the ψ profile +mat_flag = true # Construct coefficient matrices for diagnostic purposes +ode_flag = true # Integrate ODEs for stability of the internal long-wavelength mode (must be true for GPEC) +vac_flag = true # Compute plasma, vacuum, and total energies for free-boundary modes + +psiedge = 0.99 # Edge dW(ψ) diagnostic scan band [psiedge, psilim]; set ≥ psilim to disable +qlow = 1.02 # Integration initiated at q determined by min(q0, qlow) +qhigh = 1e3 # Integration terminated at q limit determined by min(qa, qhigh) +sing_start = 0 # Start integration at the sing_start'th rational from the axis (psilow) + +nn_low = 1 # Smallest toroidal mode number to include +nn_high = 1 # Largest toroidal mode number to include +delta_mlow = 8 # Expands lower bound of Fourier harmonics +delta_mhigh = 8 # Expands upper bound of Fourier harmonics +mthvac = 960 # Number of points used in splines over poloidal angle at the plasma-vacuum interface + +kinetic_source = "fixed" # Kinetic matrix source: "fixed" test matrices, or "calculated" from the kinetic NTV model +kinetic_factor = 0.0 # Scaling of kinetic matrices (0 = ideal path; >0 enables kinetic mode) +eulerlagrange_tolerance = 1e-7 # Relative tolerance for ODE integration of Euler-Lagrange equations +singfac_min = 1e-4 # Fractional distance from rational q at which ideal jump enforced +ucrit = 1e3 # Column-norm threshold that triggers solution renormalization +force_wv_symmetry = true # Enforce symmetry of the vacuum energy matrix +save_interval = 3 # Save every Nth ODE step (1=all). Always saves near rational surfaces. + +# Δ' BVP + parallel integration (see ForceFreeStatesControl docstring for details) +use_parallel = true # Run parallel FM-propagator BVP path (unlocks singular/delta_prime_matrix) +parallel_threads = 2 # BVP thread cap (1 = serial/bit-deterministic; 2 ≈ +20% speedup; ≥3 saturates) +populate_dense_xi = true # Append serial-EL pass so dense ξ is stored — REQUIRED with a [PerturbedEquilibrium] section +set_psilim_via_dmlim = false # FALSE for limited/analytical equilibria — rationals sparse, dmlim would chop too much edge +dmlim = 0.2 # Truncate integration at (last_rational_q + dmlim)/n (used when set_psilim_via_dmlim = true) +# Solovev analytic equilibrium parameters (eq_type = "sol"); see SolovevConfig in src/Equilibrium. +[SOL_INPUT] +mr = 128 # Number of radial grid zones +mz = 128 # Number of axial grid zones +ma = 128 # Number of flux grid zones +e = 1.6 # Elongation +a = 0.33 # Minor radius +r0 = 1.0 # Major radius +q0 = 1.9 # Safety factor at the magnetic axis (O-point) +p0fac = 1 # Scale on-axis pressure (β changes; Φ, q fixed) +b0fac = 1 # Scale toroidal field at constant β (Bt changes; shape, β fixed) +f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, Bp fixed) + +[KineticForces] +# Multi-ion NTV: the analytic Solovev main ion is modelled as a 50/50 D-T mix (n_i = total), +# so the summed NTV runs over D + T + the neutrality-closing impurity + electrons. +kinetic_file = "kinetic.dat" +electron = true +zimp = 6 +mimp = 12 + +[[KineticForces.ion_species]] +z = 1 +m = 2 +fraction = 0.5 + +[[KineticForces.ion_species]] +z = 1 +m = 3 +fraction = 0.5 diff --git a/examples/Solovev_kinetic_multiion_example/kinetic.dat b/examples/Solovev_kinetic_multiion_example/kinetic.dat new file mode 100644 index 00000000..96fd9fa3 --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/kinetic.dat @@ -0,0 +1,23 @@ +# Kinetic profiles consistent with this Solovev equilibrium's pressure. +# Constraint: P(psi_n) = n_e*T_e + n_i*T_i must equal the equilibrium pressure +# P(psi_n) = P0*(1 - psi_n), P0 = 4.273e4 Pa (linear, zero at the edge). +# Construction (quasineutral, n_e=n_i, T_e=T_i): +# T(psi_n) = 100 + 2600*(1 - psi_n) eV (core 2700 eV, edge 100 eV) +# n(psi_n) = P0*(1 - psi_n) / (2*e*T) (core ~4.9e19; product n0*T0 = P0/2e fixed) +# Density carries the edge falloff to a small floor at psi_n=1 so T stays finite +# (avoids T->0 in the collisionality / drift-frequency formulas); the residual edge +# pressure from the floor is ~1e-6 of P0. omega_E is a free (non-pressure) profile. +psi_n n_i_m3 n_e_m3 T_i_eV T_e_eV omega_E_rads +0.00 4.939e19 4.939e19 2700.0 2700.0 1.0e4 +0.10 4.919e19 4.919e19 2440.0 2440.0 9.0e3 +0.20 4.894e19 4.894e19 2180.0 2180.0 8.0e3 +0.30 4.862e19 4.862e19 1920.0 1920.0 7.0e3 +0.40 4.820e19 4.820e19 1660.0 1660.0 6.0e3 +0.50 4.763e19 4.763e19 1400.0 1400.0 5.0e3 +0.60 4.679e19 4.679e19 1140.0 1140.0 4.0e3 +0.70 4.546e19 4.546e19 880.0 880.0 3.0e3 +0.80 4.302e19 4.302e19 620.0 620.0 2.0e3 +0.90 3.705e19 3.705e19 360.0 360.0 1.0e3 +0.95 2.899e19 2.899e19 230.0 230.0 5.0e2 +0.99 1.058e19 1.058e19 126.0 126.0 1.0e2 +1.00 2.000e18 2.000e18 100.0 100.0 0.0 diff --git a/examples/Solovev_kinetic_multiion_example/run_example.jl b/examples/Solovev_kinetic_multiion_example/run_example.jl new file mode 100644 index 00000000..fa0a9e53 --- /dev/null +++ b/examples/Solovev_kinetic_multiion_example/run_example.jl @@ -0,0 +1,4 @@ +using Pkg; +Pkg.activate(joinpath(@__DIR__, "../..")) +using GeneralizedPerturbedEquilibrium +GeneralizedPerturbedEquilibrium.main([dirname(@__FILE__)]) diff --git a/regression-harness/cases/solovev_kinetic_multiion.toml b/regression-harness/cases/solovev_kinetic_multiion.toml new file mode 100644 index 00000000..df97e674 --- /dev/null +++ b/regression-harness/cases/solovev_kinetic_multiion.toml @@ -0,0 +1,67 @@ +# Regression case: Solovev n=1 + PE + MULTI-ION NTV torque. The Solovev main ion is modelled as a +# 50/50 D-T mix, so the KineticForces quadrature sums over D + T + electrons (this fixture is a pure +# plasma, n_e = n_i, so no impurity species is present). Tracks the summed total AND each per-species +# contribution, pinning the multi-species loop (resolve_ntv_species + combine_species_states) and the +# full-composition Zeff. The impurity-species path is covered by the unit tests + the TC-24 case. +[case] +name = "solovev_kinetic_multiion" +description = "Solovev n=1, ideal + PE + multi-ion NTV (D + T + electron summed)" +example_dir = "examples/Solovev_kinetic_multiion_example" + +# Summed total NTV over all species — the headline multi-ion quantity. +[quantities.ntv_total] +h5path = "kinetic_forces/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV total (D+T+imp+e) [N·m]" +noise_threshold = 1e-8 +order = 10 + +[quantities.ntv_D] +h5path = "kinetic_forces_ion1_z1_m2/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV Deuterium [N·m]" +noise_threshold = 1e-8 +order = 11 + +[quantities.ntv_T] +h5path = "kinetic_forces_ion2_z1_m3/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV Tritium [N·m]" +noise_threshold = 1e-8 +order = 12 + +[quantities.ntv_electron] +h5path = "kinetic_forces_electron/fgar/total_torque" +type = "real_scalar" +extract = "value" +label = "NTV electron [N·m]" +noise_threshold = 1e-8 +order = 14 + +# Stability anchors — confirm the upstream FFS/PE stages feeding the NTV diagnostic. +[quantities.et_real] +h5path = "FreeBoundaryStability/eigenmode_energies" +type = "complex_vector" +extract = "real_first" +label = "total energy Re(et[1])" +noise_threshold = 1e-10 +order = 20 + +[quantities.msing] +h5path = "singular/msing" +type = "int_scalar" +extract = "value" +label = "# singular surfaces" +noise_threshold = 0 +order = 21 + +[quantities.runtime] +h5path = "" +type = "runtime" +extract = "value" +label = "Runtime (s)" +noise_threshold = 0.0 +order = 90 diff --git a/src/Equilibrium/KineticProfiles.jl b/src/Equilibrium/KineticProfiles.jl index 606e2449..359d26c0 100644 --- a/src/Equilibrium/KineticProfiles.jl +++ b/src/Equilibrium/KineticProfiles.jl @@ -12,6 +12,42 @@ shimming. using DelimitedFiles using HDF5 +# Physical constants and collisionality normalizations shared by `load_kinetic_profiles` and +# `resolve_ntv_species` (Krook collision frequency, Logan & Park 2013 Eq. 6). +const _EV_J = 1.602e-19 # eV → J +const _KEV_J = 1.602e-16 # 1 keV in J (temperature normalization in the collision frequency) +const _NU_PREFAC = 3.5e17 # Krook collision-frequency prefactor +const _MP = 1.672_614e-27 # proton mass [kg] +const _ME = 9.109_1e-31 # electron mass [kg] + +# NRL-formulary Coulomb logarithm (natural log); n_e in m⁻³, T_e in J. +_coulomb_log(ne, Te) = 17.3 - 0.5 * log(ne / 1.0e20) + 1.5 * log(Te / _KEV_J) + +# Hard stop on an unphysical negative density (bad input or a cubic-resample overshoot). +function _assert_nonneg_density(kinetic_file, label, arr, psi_reg) + i = findfirst(<(0), arr) + i === nothing || + error("kinetic file '$kinetic_file': negative $label = $(arr[i]) at ψ_n = $(psi_reg[i]) — unphysical") +end + +""" + ResolvedNTVSpecies{P} + +One fully-resolved species in the NTV sum — a main ion, the neutrality-closing impurity, or the +electrons. Fields: charge `z`, mass `m` (proton masses), `electron` flag, `label`, and `profiles` +(a `KineticProfileSplines` view carrying this species' density as `ni_spline` and its +full-composition collision frequency as `nui_spline`). Produced by `resolve_ntv_species`; both NTV +paths — the KineticForces ψ-quadrature and the self-consistent kinetic-matrix build — iterate this +and sum τ = Σ_s τ_s. +""" +struct ResolvedNTVSpecies{P} + z::Int + m::Int + electron::Bool + label::String + profiles::P +end + """ KineticProfileData @@ -21,15 +57,17 @@ present; every profile field is `nothing` when the source omits it, so each consumer validates only the fields it needs. Units: densities m⁻³, temperatures eV, frequencies rad/s, diffusivities m²/s. -| field | meaning | -|-------------|--------------------------------------| -| `n_i` | main-ion density | -| `n_e` | electron density | -| `T_i`/`T_e` | ion / electron temperature | -| `omega_E` | ExB rotation | -| `omega_tor` | toroidal rotation (optional) | -| `chi_e` | perpendicular heat diffusivity χ⊥ | -| `chi_phi` | toroidal momentum diffusivity χ_φ | +| field | meaning | +|:------------------- |:----------------------------------------------------------------------------------------------------------------- | +| `n_i` | main-ion density | +| `n_e` | electron density | +| `T_i`/`T_e` | ion / electron temperature | +| `omega_E` | ExB rotation | +| `omega_tor` | toroidal rotation (optional) | +| `chi_e` | perpendicular heat diffusivity χ⊥ | +| `chi_phi` | toroidal momentum diffusivity χ_φ | +| `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 | """ struct KineticProfileData psi::Vector{Float64} @@ -41,14 +79,16 @@ struct KineticProfileData omega_tor::Union{Nothing,Vector{Float64}} chi_e::Union{Nothing,Vector{Float64}} chi_phi::Union{Nothing,Vector{Float64}} + species_densities::Union{Nothing,Dict{String,Vector{Float64}}} provenance::String end 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_phi=nothing, provenance="") + omega_E=nothing, omega_tor=nothing, chi_e=nothing, chi_phi=nothing, + species_densities=nothing, provenance="") _v(x) = x === nothing ? nothing : Float64.(collect(x)) 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_phi), String(provenance)) + _v(omega_E), _v(omega_tor), _v(chi_e), _v(chi_phi), species_densities, String(provenance)) end const _KINETIC_H5_EXTS = (".h5", ".hdf5", ".he5") @@ -72,15 +112,19 @@ function read_kinetic_file(path::AbstractString; group::AbstractString="/") return ext in _KINETIC_H5_EXTS ? _read_kinetic_h5(path; group=group) : _read_kinetic_ascii(path) end -"""Read a legacy 6-column ASCII kinetic table into a `KineticProfileData`.""" +""" +Read a legacy 6-column ASCII kinetic table 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))") end -"""Read the GPEC HDF5 kinetic schema. Only `psi` is required; other datasets are -optional and surfaced as-is.""" +""" +Read the GPEC HDF5 kinetic schema. Only `psi` is required; other datasets are +optional and surfaced as-is. +""" function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") h5open(path, "r") do f g = group == "/" ? f : f[group] @@ -89,10 +133,20 @@ function _read_kinetic_h5(path::AbstractString; group::AbstractString="/") prov = "" ats = attributes(g) haskey(ats, "provenance") && (prov = string(read(ats["provenance"]))) + # Any dataset outside the standard schema is treated as a named per-species density + # profile (e.g. "n_D", "n_T") for explicit multi-ion input. + standard = ("psi", "n_i", "n_e", "T_i", "T_e", "omega_E", "omega_tor", "chi_e", "chi_phi") + extras = Dict{String,Vector{Float64}}() + for k in keys(g) + k in standard && continue + v = read(g[k]) + v isa AbstractArray && (extras[k] = Float64.(vec(v))) + end return KineticProfileData(; psi=Float64.(vec(read(g["psi"]))), n_i=rd("n_i"), n_e=rd("n_e"), T_i=rd("T_i"), T_e=rd("T_e"), omega_E=rd("omega_E"), omega_tor=rd("omega_tor"), - chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), provenance=prov) + chi_e=rd("chi_e"), chi_phi=rd("chi_phi"), + species_densities=(isempty(extras) ? nothing : extras), provenance=prov) end end @@ -128,6 +182,145 @@ function write_kinetic_h5(path::AbstractString, data::KineticProfileData; return path end +""" + multi_ion_composition(zs, ns, ne; zimp, mimp) -> (zeff, zpitch, n_main, n_imp) + +Full-composition effective charge and pitch-angle enhancement for a plasma with an +arbitrary number of main-ion species. Inputs are per-point (scalars or broadcast arrays): + + - `zs`, `ns`: main-ion charges and number densities (one entry per species) + - `ne`: electron density + - `zimp`, `mimp`: single trailing impurity charge / mass; its density closes quasineutrality + +Composition (per point): + + n_main = Σ_s n_s total main-ion number density + n_imp = (n_e − Σ_s z_s n_s) / z_imp impurity density from quasineutrality + Zeff = (Σ_s z_s² n_s + z_imp² n_imp) / n_e + zpitch = 1 + (1+m_imp)/(2 m_imp)·z_imp·(Zeff−1)/(z_imp−Zeff) momentum-restoring closure + +Reduces **exactly** to the single-ion form `Zeff = z_imp − (n_i/n_e)·z_i·(z_imp−z_i)` for one +z-species (verified algebraically). The per-species pitch-angle collision frequency built from +this is `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s · (T_s)^{3/2})` — shared `zpitch`, `n_main`, +`Zeff`, `lnΛ`; the test species contributes only its own `z_s²`, `m_s`, `T_s`. +""" +function multi_ion_composition(zs::AbstractVector, ns::AbstractVector, ne::Real; + zimp::Real, mimp::Real) + n_main = sum(ns) + charge_density = sum(z * n for (z, n) in zip(zs, ns)) # Σ z_s n_s + n_imp = ne > 0 ? (ne - charge_density) / zimp : 0.0 + zeff = ne > 0 ? (sum(z^2 * n for (z, n) in zip(zs, ns)) + zimp^2 * n_imp) / ne : Float64(zimp) + zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (zeff - 1.0) / (zimp - zeff) + return (zeff, zpitch, n_main, n_imp) +end + +""" + resolve_ntv_species(kinetic_file, ion_species; electron, zimp, mimp, ...) -> Vector{ResolvedNTVSpecies} + +Resolve the **full NTV species set** for a multi-ion run: the main ions, the neutrality-closing +impurity, and (if `electron`) the electrons. Returns `ResolvedNTVSpecies` descriptors — the +single set that BOTH NTV paths (the KineticForces quadrature and the self-consistent kinetic-matrix +build) loop over and sum. Errors on any negative density (unphysical input / resample overshoot). + +Each descriptor's `profiles` view carries **that species' resonant density** `n_s` (as `ni_spline`) +and its **full-composition collision frequency** `ν_s` (as `nui_spline`); `ne/Te/ωE/Zeff/lnΛ` (and +the electron `nue`) are shared. Main-ion density is resolved from `IonSpecies`: `fraction` ⇒ +`n_s = fraction · n_i` (the kinetic file's `n_i` = total main-ion density); explicit per-species +profiles (`density`) are not yet wired. Composition (Zeff, zpitch, `n_imp`) comes from +`multi_ion_composition`; per-species `ν_s = (zpitch/3.5e17)·z_s²·n_main·lnΛ / (√m_s·(T_i)^{3/2})`. +The impurity (`zimp`,`mimp`, density = the quasineutrality `n_imp`) is included as its own resonant +species. Reduces to `load_kinetic_profiles` for one z=1 main ion with `fraction=1`. +""" +function resolve_ntv_species(kinetic_file::AbstractString, ion_species::AbstractVector; + electron::Bool=false, zimp::Integer=6, mimp::Integer=12, + density_factor::Float64=1.0, temperature_factor::Float64=1.0, + ExB_rotation_factor::Float64=1.0, toroidal_rotation_factor::Float64=1.0) + + (density_factor == 1.0 && temperature_factor == 1.0 && ExB_rotation_factor == 1.0 && + toroidal_rotation_factor == 1.0) || + error("profile scaling factors are not yet supported together with a multi-ion ion_species list") + + data = read_kinetic_file(kinetic_file) + _need(f, n) = f === nothing ? error("kinetic file '$kinetic_file' missing required '$n'") : f + psi_in = data.psi + ne_in = _need(data.n_e, "n_e") + Ti_in = _need(data.T_i, "T_i") + Te_in = _need(data.T_e, "T_e") + omE_in = _need(data.omega_E, "omega_E") + ni_in = data.n_i === nothing ? copy(ne_in) : data.n_i # TOTAL main-ion density + + nkin = 100 + psi_reg = collect(0:nkin) ./ nkin + ne = _cubic_resample(psi_in, ne_in, psi_reg) + ni_total = _cubic_resample(psi_in, ni_in, psi_reg) + Ti = _cubic_resample(psi_in, Ti_in, psi_reg) .* _EV_J + Te = _cubic_resample(psi_in, Te_in, psi_reg) .* _EV_J + omegaE = _cubic_resample(psi_in, omE_in, psi_reg) + _assert_nonneg_density(kinetic_file, "n_e", ne, psi_reg) + _assert_nonneg_density(kinetic_file, "n_i (total)", ni_total, psi_reg) + + # Resolve each species' resonant density: `fraction` ⇒ share of the total n_i; `density` ⇒ + # an explicit named profile from the (HDF5) kinetic file (resampled to the working grid). + zs = [Int(s.z) for s in ion_species] + ns = Vector{Vector{Float64}}(undef, length(ion_species)) + for (si, s) in enumerate(ion_species) + has_frac = !isnan(s.fraction) + has_dens = !isempty(s.density) + (has_frac ⊻ has_dens) || error("ion_species[$si]: specify exactly one of `fraction` or `density`") + if has_dens + (data.species_densities !== nothing && haskey(data.species_densities, s.density)) || + error( + "ion_species[$si]: density profile `$(s.density)` not found in kinetic file " * + "(available: $(data.species_densities === nothing ? "none — ASCII files carry only n_i" : join(keys(data.species_densities), ", ")))" + ) + ns[si] = _cubic_resample(psi_in, data.species_densities[s.density], psi_reg) + else + ns[si] = s.fraction .* ni_total + end + _assert_nonneg_density(kinetic_file, "ion_species[$si] density", ns[si], psi_reg) + end + + # Shared composition + collisionality (natural-log Coulomb log), per grid point. + npts = length(psi_reg) + zeff = zeros(npts) + zpitch = zeros(npts) + loglam = zeros(npts) + n_main = zeros(npts) + n_imp = zeros(npts) + nue = zeros(npts) + 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] = (ne[i] > 0 && Te[i] > 0) ? _coulomb_log(ne[i], Te[i]) : 0.0 + nue[i] = Te[i] > 0 ? (zp / _NU_PREFAC) * ne[i] * loglam[i] / (sqrt(_ME / _MP) * (Te[i] / _KEV_J)^1.5) : 0.0 + end + # A negative impurity density means the main ions over-neutralize (Σ z_s n_s > n_e) — unphysical. + _assert_nonneg_density(kinetic_file, "impurity density (n_e < Σ z_s n_s)", n_imp, psi_reg) + + # Per-species ν_s: shared zpitch/n_main/Zeff/lnΛ, species-specific z²/m/T_i (Krook deflection + # frequency, test-particle z², Logan & Park 2013 Eq. 6). `zpitch` is a main-ion closure, applied + # approximately to the impurity/electron test species. + _nu(zsp, msp) = [Ti[i] > 0 ? (zpitch[i] / _NU_PREFAC) * zsp^2 * n_main[i] * loglam[i] / (sqrt(Float64(msp)) * (Ti[i] / _KEV_J)^1.5) : 0.0 for i in 1:npts] + _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. + species = ResolvedNTVSpecies[] + for (si, s) in enumerate(ion_species) + push!(species, ResolvedNTVSpecies(Int(s.z), Int(s.m), false, "ion$(si)_z$(s.z)_m$(s.m)", _view(ns[si], _nu(s.z, s.m)))) + end + if any(>(0), n_imp) + push!(species, ResolvedNTVSpecies(Int(zimp), Int(mimp), false, "impurity_z$(zimp)_m$(mimp)", _view(n_imp, _nu(zimp, mimp)))) + end + if electron + # The electron view's `ni_spline` is unused (the electron path reads `ne_spline`); carry ne. + push!(species, ResolvedNTVSpecies(-1, 1, true, "electron", _view(ne, nue))) + end + return identity.(species) # narrow to a concrete-eltype Vector{ResolvedNTVSpecies{...}} +end + """ load_kinetic_profiles(kinetic_file::AbstractString; zi::Int=1, zimp::Int=6, mi::Int=2, mimp::Int=12, @@ -162,11 +355,12 @@ six-column whitespace-separated ASCII table (header rows are filtered out): # Scaling sequence When any of `density_factor`, `temperature_factor`, `toroidal_rotation_factor` differ from 1.0: -1. Build first-pass cubic splines from unscaled profiles (for derivatives) -2. Compute diamagnetic frequencies `wdian`, `wdiat` and total toroidal rotation `wphi` -3. Scale: `wdian_new = temperature_factor * wdian`, `wdiat_new = temperature_factor * wdiat` (`density_factor` cancels in `T*(dn/dψ)/n`) -4. Reform: `omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new` -5. Scale density/temperature arrays: `ni *= density_factor`, `Ti *= temperature_factor`, etc. + + 1. Build first-pass cubic splines from unscaled profiles (for derivatives) + 2. Compute diamagnetic frequencies `wdian`, `wdiat` and total toroidal rotation `wphi` + 3. Scale: `wdian_new = temperature_factor * wdian`, `wdiat_new = temperature_factor * wdiat` (`density_factor` cancels in `T*(dn/dψ)/n`) + 4. Reform: `omegaE = toroidal_rotation_factor * wphi - wdian_new - wdiat_new` + 5. Scale density/temperature arrays: `ni *= density_factor`, `Ti *= temperature_factor`, etc. Then `ExB_rotation_factor` is applied independently: `omegaE *= ExB_rotation_factor`. @@ -192,10 +386,6 @@ function load_kinetic_profiles(kinetic_file::AbstractString; omegaE_input = _need(data.omega_E, "omega_E") ni_input = data.n_i === nothing ? copy(ne_input) : data.n_i - eV_to_J = 1.602e-19 - mp = 1.672_614e-27 - me = 9.109_1e-31 - nkin = 100 psi_reg = collect(0:nkin) ./ nkin @@ -207,8 +397,8 @@ function load_kinetic_profiles(kinetic_file::AbstractString; # denominator). See feedback_kf_kin_profile_linear_interp.md. ni = _cubic_resample(psi_input, ni_input, psi_reg) ne = _cubic_resample(psi_input, ne_input, psi_reg) - Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* eV_to_J - Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* eV_to_J + Ti = _cubic_resample(psi_input, Ti_input_eV, psi_reg) .* _EV_J + Te = _cubic_resample(psi_input, Te_input_eV, psi_reg) .* _EV_J omegaE = _cubic_resample(psi_input, omegaE_input, psi_reg) needs_rotation_reform = density_factor != 1.0 || temperature_factor != 1.0 || toroidal_rotation_factor != 1.0 @@ -228,8 +418,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString; dTi_dpsi = deriv1(Ti_spl) # Compute original wdian, wdiat, wphi and reform omegaE at each grid point - echarge = 1.602e-19 - chrg_ion = zi * echarge + chrg_ion = zi * _EV_J for i in eachindex(omegaE) ψ = psi_reg[i] wdian_i = ni[i] > 0 ? -2π * Ti[i] * dni_dpsi(ψ) / (chrg_ion * chi1 * ni[i]) : 0.0 @@ -258,7 +447,7 @@ function load_kinetic_profiles(kinetic_file::AbstractString; nue = zeros(Float64, nkin + 1) zeff = zeros(Float64, nkin + 1) - for i in 1:(nkin + 1) + for i in 1:(nkin+1) n_i = ni[i] n_e = ne[i] T_i = Ti[i] @@ -267,13 +456,15 @@ function load_kinetic_profiles(kinetic_file::AbstractString; z = n_e > 0 ? zimp - (n_i / n_e) * zi * (zimp - zi) : Float64(zimp) zpitch = 1.0 + (1.0 + mimp) / (2.0 * mimp) * zimp * (z - 1.0) / (zimp - z) - ll = 17.3 - 0.5 * log10(n_e / 1.0e20) + 1.5 * log10(T_e / 1.602e-16) + ll = _coulomb_log(n_e, T_e) loglam[i] = ll + # Test-particle z² pitch-angle (Krook deflection) scaling — Logan & Park 2013 Eq. 6; + # generalizes the implicit zi=1 form; identical to the single-species limit of resolve_ntv_species. nui[i] = T_i > 0 ? - (zpitch / 3.5e17) * n_i * ll / (sqrt(1.0 * mi) * (T_i / 1.602e-16)^1.5) : 0.0 + (zpitch / _NU_PREFAC) * zi^2 * n_i * ll / (sqrt(1.0 * mi) * (T_i / _KEV_J)^1.5) : 0.0 nue[i] = T_e > 0 ? - (zpitch / 3.5e17) * n_e * ll / (sqrt(me / mp) * (T_e / 1.602e-16)^1.5) : 0.0 + (zpitch / _NU_PREFAC) * n_e * ll / (sqrt(_ME / _MP) * (T_e / _KEV_J)^1.5) : 0.0 zeff[i] = z end diff --git a/src/ForceFreeStates/EulerLagrange.jl b/src/ForceFreeStates/EulerLagrange.jl index 9f93e60a..1f8fe3bc 100644 --- a/src/ForceFreeStates/EulerLagrange.jl +++ b/src/ForceFreeStates/EulerLagrange.jl @@ -676,7 +676,7 @@ function cross_kinetic_singular_surf!( ksurf = intr.kinsing[ising] dpsi = ksurf.psifac - odet.psifac - params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, 0)) + params = (ctrl, equil, ffit, intr, odet, IntegrationChunk(0.0, 0.0, false, ising, 1)) du1 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) du2 = zeros(ComplexF64, intr.numpert_total, intr.numpert_total, 2) diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index edee8303..104a00e0 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -218,6 +218,10 @@ function main_from_inputs( KineticForces.KineticForcesControl() kinetic_profiles = nothing + # Multi-species NTV set (main ions + neutrality impurity + electrons), resolved once and + # summed by BOTH the self-consistent kinetic-matrix build and the KineticForces quadrature. + # `nothing` ⇒ single-species from kf_ctrl.zi/mi (unchanged behaviour). + kf_species = nothing needs_kinetic_profiles = haskey(inputs, "KineticForces") || (ctrl.kinetic_factor > 0 && ctrl.kinetic_source == "calculated") if needs_kinetic_profiles @@ -229,6 +233,12 @@ function main_from_inputs( density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor, chi1=2π * equil.psio) + if !isempty(kf_ctrl.ion_species) + kf_species = Equilibrium.resolve_ntv_species(kinetic_file, kf_ctrl.ion_species; + electron=kf_ctrl.electron, zimp=kf_ctrl.zimp, mimp=kf_ctrl.mimp, + density_factor=kf_ctrl.density_factor, temperature_factor=kf_ctrl.temperature_factor, + ExB_rotation_factor=kf_ctrl.ExB_rotation_factor, toroidal_rotation_factor=kf_ctrl.toroidal_rotation_factor) + end end # Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, @@ -418,7 +428,7 @@ function main_from_inputs( calculated_cb = (c, e, i, m, f) -> KineticForces.compute_calculated_kinetic_matrices( c, e, i, m, f; - kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles) + kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=kf_species) make_kinetic_matrix(ctrl, equil, ffit, intr, metric; calculated_source=calculated_cb) @@ -605,12 +615,37 @@ function main_from_inputs( kf_intr = KineticForces.KineticForcesInternal(equil; verbose=kf_ctrl.verbose) KineticForces.set_perturbation_data!(kf_intr, pe_state, intr, equil, metric) - kf_state = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles) - - if kf_ctrl.write_outputs_to_HDF5 - h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file - KineticForces.write_to_hdf5!(h5file, kf_state) + if isempty(kf_ctrl.ion_species) + # Single main ion (from zi/mi) — unchanged path. + kf_state = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, equil, kinetic_profiles) + if kf_ctrl.write_outputs_to_HDF5 + h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file + KineticForces.write_to_hdf5!(h5file, kf_state) + end + end + else + # Multi-species: summed over main ions + neutrality impurity + electrons + # (kf_species, resolved once above and shared with the kinetic-matrix build). + states = KineticForces.KineticForcesState[] + labels = String[] + for sp in kf_species + sctrl = deepcopy(kf_ctrl) + sctrl.zi = sp.z; sctrl.mi = sp.m; sctrl.electron = sp.electron + sctrl.ion_species = KineticForces.IonSpecies[] + st = KineticForces.KineticForcesState() + KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, equil, sp.profiles) + push!(states, st); push!(labels, sp.label) + end + @info "Multi-species NTV: summed $(length(kf_species)) species ($(join((s.label for s in kf_species), ", ")))" + kf_state = KineticForces.combine_species_states(states) + if kf_ctrl.write_outputs_to_HDF5 + h5open(joinpath(intr.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file + for (st, lbl) in zip(states, labels) + KineticForces.write_to_hdf5!(h5file, st; group_name="kinetic_forces_$lbl") + end + KineticForces.write_to_hdf5!(h5file, kf_state) # "kinetic_forces" = total + end end end end diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index d36fead9..71f58f0b 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -31,22 +31,25 @@ tracked as follow-up work blocked on PR #196 — see the plan's "Out of scope" section. # Arguments -- `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`) -- `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines -- `ffs_intr`: ForceFreeStatesInternal (mode indexing) -- `metric`: MetricData (provides ψ grid via `metric.xs`) -- `ffit`: FourFitVars (used only for `numpert_total` cross-check) + + - `ffs_ctrl`: ForceFreeStatesControl (carries `kinetic_factor`, `kinetic_source`) + - `equil`: PlasmaEquilibrium with 2D interpolants and named profile/geometry splines + - `ffs_intr`: ForceFreeStatesInternal (mode indexing) + - `metric`: MetricData (provides ψ grid via `metric.xs`) + - `ffit`: FourFitVars (used only for `numpert_total` cross-check) # Keyword arguments -- `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to - carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the - KineticForces kernel needs but ForceFreeStatesControl does not expose. -- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- - profile splines loaded via `Equilibrium.load_kinetic_profiles`. + + - `kf_ctrl`: KineticForcesControl, defaults to `KineticForcesControl()`. Used to + carry NTV-specific knobs (nl, zi, mi, wdfac, divxfac, electron) that the + KineticForces kernel needs but ForceFreeStatesControl does not expose. + - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Required. Named kinetic- + profile splines loaded via `Equilibrium.load_kinetic_profiles`. # Returns -- `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)` -- `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)` + + - `kw_flat::Array{ComplexF64,3}`: Energy matrices, shape `(mpsi, np^2, 6)` + - `kt_flat::Array{ComplexF64,3}`: Torque matrices, shape `(mpsi, np^2, 6)` """ function compute_calculated_kinetic_matrices( _ffs_ctrl, @@ -54,8 +57,9 @@ function compute_calculated_kinetic_matrices( ffs_intr, metric, ffit; - kf_ctrl::KineticForcesControl = KineticForcesControl(), + kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, + species::Union{Nothing,AbstractVector}=nothing ) xs = metric.xs mpsi = length(xs) @@ -95,44 +99,53 @@ function compute_calculated_kinetic_matrices( # are read-only and safely shared through deepcopy semantics. nl = kf_ctrl.nl nthreads = Threads.maxthreadid() - thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads] - thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] - thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] + thread_intrs = [deepcopy(kf_intr) for _ in 1:nthreads] + thread_full_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] + thread_full_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] thread_block_w = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] thread_block_t = [zeros(ComplexF64, mpert, mpert, 6) for _ in 1:nthreads] - Threads.@threads for ipsi in 1:mpsi - tid = Threads.threadid() - intr_t = thread_intrs[tid] - full_w = thread_full_w[tid] - full_t = thread_full_t[tid] - block_w = thread_block_w[tid] - block_t = thread_block_t[tid] - psi = xs[ipsi] - for in_idx in 1:npert - n = ffs_intr.nlow + in_idx - 1 - fill!(full_w, 0) - fill!(full_t, 0) - for ell in -nl:nl - fill!(block_w, 0) - fill!(block_t, 0) - compute_kinetic_matrices_at_psi!( - block_w, block_t, psi, n, ell, - kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.wdfac, kf_ctrl.divxfac, - kf_ctrl.electron, equil, intr_t, kinetic_profiles; - nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, - atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda - ) - full_w .+= block_w - full_t .+= block_t - end - - # Place the n-block on the diagonal of the full np×np matrix and flatten. - row_offset = (in_idx - 1) * mpert - for k in 1:6, j in 1:mpert, i in 1:mpert - idx = (row_offset + j - 1) * np + (row_offset + i) - kw_flat[ipsi, idx, k] = full_w[i, j, k] - kt_flat[ipsi, idx, k] = full_t[i, j, k] + # Species set summed into the kinetic matrices: the resolved multi-species set + # (main ions + neutrality impurity + electrons) when provided, else a single species + # from kf_ctrl. The kinetic W/torque matrices are additive over species, so we + # accumulate (+=) each species' block-diagonal contribution. + splist = species === nothing ? + [Equilibrium.ResolvedNTVSpecies(kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.electron, "single", kinetic_profiles)] : species + + for sp in splist + Threads.@threads for ipsi in 1:mpsi + tid = Threads.threadid() + intr_t = thread_intrs[tid] + full_w = thread_full_w[tid] + full_t = thread_full_t[tid] + block_w = thread_block_w[tid] + block_t = thread_block_t[tid] + psi = xs[ipsi] + for in_idx in 1:npert + n = ffs_intr.nlow + in_idx - 1 + fill!(full_w, 0) + fill!(full_t, 0) + for ell in -nl:nl + fill!(block_w, 0) + fill!(block_t, 0) + compute_kinetic_matrices_at_psi!( + block_w, block_t, psi, n, ell, + sp.z, sp.m, kf_ctrl.wdfac, kf_ctrl.divxfac, + sp.electron, equil, intr_t, sp.profiles; + nutype=kf_ctrl.nutype, f0type=kf_ctrl.f0type, nufac=kf_ctrl.nufac, + atol_xlmda=kf_ctrl.atol_xlmda, rtol_xlmda=kf_ctrl.rtol_xlmda + ) + full_w .+= block_w + full_t .+= block_t + end + + # Place the n-block on the diagonal of the full np×np matrix and accumulate. + row_offset = (in_idx - 1) * mpert + for k in 1:6, j in 1:mpert, i in 1:mpert + idx = (row_offset + j - 1) * np + (row_offset + i) + kw_flat[ipsi, idx, k] += full_w[i, j, k] + kt_flat[ipsi, idx, k] += full_t[i, j, k] + end end end end diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 4175a62b..d2ce59ec 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -41,7 +41,7 @@ Warn when the ψ torque quadrature terminated without satisfying the requested t silent-garbage scenario for weak applied fields, since NTV scales as δB². """ function check_psi_quadrature_convergence(total::ComplexF64, quad_err::Float64, - ctrl::KineticForcesControl, method::String) + ctrl::KineticForcesControl, method::String) if quad_err > max(ctrl.atol_psi, ctrl.rtol_psi * abs(total)) @warn "ψ torque quadrature ($method) did not converge within maxevals_psi=$(ctrl.maxevals_psi): " * "error estimate $quad_err vs |T|=$(abs(total)) N·m. Raise maxevals_psi or loosen rtol_psi." @@ -63,14 +63,16 @@ diagnostic T(ψ) profile at no extra cost (the values are computed anyway — we just keep them). # Returns + NamedTuple with: -- `total::ComplexF64`: Total integrated torque -- `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points -- `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method) -- `psi_nsteps::Int`: Number of integrand evaluations -- `psi_quad_error::Float64`: Quadrature error estimate for the total torque -- `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces) -- `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting + + - `total::ComplexF64`: Total integrated torque + - `torque_profile`: NamedTuple of (psi, dtdpsi, t_cumulative) from evaluation points + - `matrix_integrated`: Trapezoidal-integrated mpert×mpert×6 matrix (if matrix method) + - `psi_nsteps::Int`: Number of integrand evaluations + - `psi_quad_error::Float64`: Quadrature error estimate for the total torque + - `panel_psis::Vector{Float64}`: Quadrature panel boundaries actually used (bounds + interior resonant surfaces) + - `resonance_psis::Vector{Float64}`: Located kinetic-resonance ψ surfaces (Ω_ℓ(x=1)=0), for diagnostics/plotting The integral is paneled at the rational-surface ψ locations (`intr.sing_psis`) and capped at `ctrl.maxevals_psi` evaluations; a warning is emitted if the quadrature fails to reach @@ -95,7 +97,7 @@ function integrate_psi_quadgk( if x0 >= xout return (total=ComplexF64(0.0), torque_profile=nothing, matrix_integrated=nothing, psi_nsteps=0, psi_quad_error=0.0, - panel_psis=Float64[], resonance_psis=Float64[]) + panel_psis=Float64[], resonance_psis=Float64[]) end # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine @@ -133,9 +135,9 @@ function integrate_psi_quadgk( w = is_matrix_method ? thread_wtw[tid] : nothing is_matrix_method && fill!(w, 0) tpsi!(thread_tpsi[tid], psi, n, l, zi, mi, wdfac, divxfac, - electron, method, equil, thread_intrs[tid], kinetic_profiles; - op_wmats=w, - atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) + electron, method, equil, thread_intrs[tid], kinetic_profiles; + op_wmats=w, + atol_xlmda=ctrl.atol_xlmda, rtol_xlmda=ctrl.rtol_xlmda) harm_vals[ell_idx] = thread_tpsi[tid][] is_matrix_method && (harm_elems[ell_idx] .= w) end @@ -205,7 +207,7 @@ function integrate_psi_quadgk( "$(length(pts) - 1) panels ($n_rational rational + $n_resonance kinetic resonance surfaces)" return (total=total, torque_profile=torque_profile, matrix_integrated=matrix_integrated, psi_nsteps=npts, psi_quad_error=quad_err, - panel_psis=pts, resonance_psis=sort(resonance_psis)) + panel_psis=pts, resonance_psis=sort(resonance_psis)) end @@ -213,6 +215,61 @@ end # High-level orchestration # ============================================================================ +""" + combine_species_states(states) -> KineticForcesState + +Sum per-species [`KineticForcesState`](@ref) results into a single total (τ = Σ_s τ_s). Per +method: the scalar `total_torque`/`total_energy` are summed exactly; the dT/dψ profile is summed +by linearly interpolating each species' own `(psi_grid, dtdpsi)` arrays onto the sorted union of +the species ψ grids (zero outside a species' range), and the cumulative T(ψ) is re-integrated +(trapezoid) from that summed profile. All species share the same ψ-integration range +(`ctrl.psilims`), so the grids differ only in adaptive nodes. The combined `MethodResult` carries +only the summed scalars and profile; per-species diagnostics (`torque_profile`, `records`, +`panel_psis`, `resonance_psis`) are not aggregated and are left at their defaults. +""" +function combine_species_states(states::AbstractVector{KineticForcesState}) + combined = KineticForcesState() + isempty(states) && return combined + method_names = collect(keys(first(states).method_results)) + for mname in method_names + results = [s.method_results[mname] for s in states if haskey(s.method_results, mname)] + isempty(results) && continue + grid = sort(unique(reduce(vcat, (r.psi_grid for r in results); init=Float64[]))) + # Sum each species' dT/dψ onto the union grid by linear interpolation of its own + # (psi_grid, dtdpsi) arrays (torque_profile interpolants are not populated here); zero + # outside a species' ψ range. + dtdpsi = zeros(ComplexF64, length(grid)) + for r in results + length(r.psi_grid) >= 2 || continue + o = sortperm(r.psi_grid) + xs = r.psi_grid[o] + ys = r.dtdpsi[o] + for (k, ψ) in enumerate(grid) + (ψ < xs[1] || ψ > xs[end]) && continue + j = searchsortedlast(xs, ψ) + if j == length(xs) + dtdpsi[k] += ys[end] + else + t = (ψ - xs[j]) / (xs[j+1] - xs[j]) + dtdpsi[k] += (1 - t) * ys[j] + t * ys[j+1] + end + end + end + tcum = zeros(ComplexF64, length(grid)) + for j in 2:length(grid) + tcum[j] = tcum[j-1] + 0.5 * (dtdpsi[j] + dtdpsi[j-1]) * (grid[j] - grid[j-1]) + end + combined.method_results[mname] = MethodResult(; + method=mname, nn=first(results).nn, + total_torque=sum(r.total_torque for r in results), + total_energy=sum(r.total_energy for r in results), + psi_grid=grid, dtdpsi=dtdpsi, t_cumulative=tcum, + psi_nsteps=sum(r.psi_nsteps for r in results)) + end + combined.completed = true + return combined +end + """ compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, ctrl::KineticForcesControl, equil, kinetic_profiles) @@ -224,15 +281,16 @@ For multi-n calculations, loops over toroidal mode numbers and assembles block-diagonal kinetic matrices. # Arguments -- `state::KineticForcesState`: Accumulates results for all methods -- `intr::KineticForcesInternal`: Internal state with equilibrium data -- `ctrl::KineticForcesControl`: Control parameters specifying which methods to run -- `equil`: PlasmaEquilibrium with 2D interpolants -- `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines + + - `state::KineticForcesState`: Accumulates results for all methods + - `intr::KineticForcesInternal`: Internal state with equilibrium data + - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run + - `equil`: PlasmaEquilibrium with 2D interpolants + - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines """ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, - ctrl::KineticForcesControl, equil, - kinetic_profiles::Equilibrium.KineticProfileSplines) + ctrl::KineticForcesControl, equil, + kinetic_profiles::Equilibrium.KineticProfileSplines) for entry in METHOD_REGISTRY getfield(ctrl, entry.flag) || continue @@ -311,7 +369,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor t_cumulative=t_cum_out, psi_nsteps=psi_nsteps_total, panel_psis=panel_psis_out, - resonance_psis=resonance_psis_out, + resonance_psis=resonance_psis_out ) state.method_results[method] = result_entry diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 7e533d8e..9dcc6a58 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -3,15 +3,16 @@ Single source of truth for the NTV calculation methods. Each entry is a NamedTuple `(name, flag, kind, doc)`: -- `name` — short method identifier used as the HDF5 group key and in `intr.method` -- `flag` — the `KineticForcesControl` field symbol that enables the method -- `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` - (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the - three special-cased methods) -- `doc` — one-line description printed in verbose output + + - `name` — short method identifier used as the HDF5 group key and in `intr.method` + - `flag` — the `KineticForcesControl` field symbol that enables the method + - `kind` — dispatch routing tag consumed by `method_kind` / `Torque.jl` + (`:gar` for the GAR/matrix family, `:fcgl`/`:rlar`/`:clar` for the + three special-cased methods) + - `doc` — one-line description printed in verbose output The method names/docs and the `Compute.jl` enable list are all derived from this -tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. +tuple, and `Torque.jl` routes on `kind`, so the methods are enumerated in one place. To add a method: append an entry here and add the matching `*_flag` field to `KineticForcesControl`. """ @@ -49,13 +50,48 @@ function method_kind(name::AbstractString) end +""" + IonSpecies(; z, m, fraction=NaN, density="") + +One main-ion species in a multi-ion NTV run. `z`/`m` are the charge (e) and mass (proton +masses). The density is given by **exactly one** of: + + - `fraction`: this species' share of the total main-ion density (the kinetic file's `n_i` + column), e.g. 0.5 for each of 50/50 D-T. Convenient when species share a profile shape. + - `density`: the name of an explicit per-species density profile in the (HDF5) kinetic file, + for measured profiles with distinct shapes. + +Both resolve to a per-species density profile `n_s(ψ)` at load time (one downstream path). +Constructed from a TOML `[[KineticForces.ion_species]]` table via `convert`. +""" +struct IonSpecies + z::Int + m::Int + fraction::Float64 + density::String +end +IonSpecies(; z::Integer, m::Integer, fraction=NaN, density="") = + IonSpecies(Int(z), Int(m), Float64(fraction), String(density)) +IonSpecies(d::AbstractDict) = IonSpecies(; (Symbol(k) => v for (k, v) in d)...) +# TOML.jl parses [[KineticForces.ion_species]] as a Vector{Any} of Dicts, so accept any +# AbstractVector and build each element (pass IonSpecies through unchanged). +Base.convert(::Type{Vector{IonSpecies}}, v::AbstractVector) = + IonSpecies[x isa IonSpecies ? x : IonSpecies(x) for x in v] + """ KineticForcesControl User-facing control parameters from the TOML `[KineticForces]` section. Configures which NTV methods to run, species parameters, tolerances, and output options. +`ion_species` selects the plasma model: empty (default) ⇒ a single main ion from `zi`/`mi` +(unchanged behaviour); non-empty ⇒ a multi-main-ion run where the NTV is computed per +`IonSpecies` under one shared full-composition Zeff and summed. For a multi-ion run the kinetic +file's `n_i` column is the **total** main-ion density, split across species by their `fraction` +or explicit `density` profile. + Constructed via keyword arguments or from a TOML dict: + ```julia ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` @@ -90,6 +126,7 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce zimp::Int = 6 # Impurity charge mimp::Int = 12 # Impurity mass electron::Bool = false # Include electron contribution + ion_species::Vector{IonSpecies} = IonSpecies[] # multi-main-ion set (see docstring); empty ⇒ single ion # Mode numbers nn::Int = 1 # Toroidal mode number @@ -176,12 +213,13 @@ Internal working state for KineticForces calculations. Holds equilibrium-derived quantities, profile interpolants, and integration results. Fields replacing former module-level globals: -- `ro`, `bo`, `chi1`: Equilibrium geometry parameters -- `mthsurf`, `mfac`: Poloidal grid info -- `dbob_m`, `divx_m`: Perturbation mode interpolants -- `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as - panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on - Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection + + - `ro`, `bo`, `chi1`: Equilibrium geometry parameters + - `mthsurf`, `mfac`: Poloidal grid info + - `dbob_m`, `divx_m`: Perturbation mode interpolants + - `sing_psis`: Rational-surface ψ locations (sorted, from the stability analysis), used as + panel boundaries for the outer ψ torque quadrature so the resonant peaks fall on + Gauss-Kronrod interval endpoints instead of driving deep adaptive bisection Equilibrium and kinetic profile data are read directly from the `PlasmaEquilibrium` (`equil.profiles`, `equil.geometry`) and the @@ -271,17 +309,17 @@ function KineticForcesInternal(equil; verbose::Bool=false) # Axis toroidal field F(0)/ro that normalizes λ = μ·bo/E; F_spline stores 2πF. bo_axis = abs(equil.profiles.F_spline(0.0)) / (2π * equil.ro) KineticForcesInternal(; - ro = equil.ro, - bo = bo_axis, - chi1 = 2π * equil.psio, + ro=equil.ro, + bo=bo_axis, + chi1=2π * equil.psio, mthsurf, - tpsi_xs = collect(range(0.0, 1.0, length=nth)), - tpsi_B = Vector{Float64}(undef, nth), - tpsi_dBdpsi = Vector{Float64}(undef, nth), - tpsi_dBdtheta = Vector{Float64}(undef, nth), - tpsi_jac = Vector{Float64}(undef, nth), - tpsi_djdpsi = Vector{Float64}(undef, nth), - verbose, + tpsi_xs=collect(range(0.0, 1.0; length=nth)), + tpsi_B=Vector{Float64}(undef, nth), + tpsi_dBdpsi=Vector{Float64}(undef, nth), + tpsi_dBdtheta=Vector{Float64}(undef, nth), + tpsi_jac=Vector{Float64}(undef, nth), + tpsi_djdpsi=Vector{Float64}(undef, nth), + verbose ) end @@ -291,19 +329,21 @@ end Populate perturbation data from PerturbedEquilibriumState into KineticForcesInternal. Builds three interpolant sets from PE Clebsch displacements: -1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ -2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) -3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting + + 1. `xs_m` — [ξ^ψ, ∂ξ^ψ/∂ψ, ξ^α] CubicSeriesInterpolants over ψ + 2. `dbob_m` — δB/B Fourier modes via JBB deweighting (Fortran set_peq) + 3. `divx_m` — ∇·ξ⊥ Fourier modes via JBB deweighting The JBB deweighting algorithm (Fortran pentrc/inputs.f90:828-868): -1. Apply geometric matrices S,T,X,Y,Z in m-space -2. Inverse DFT to θ-space -3. Divide by J·B² at each θ -4. Forward DFT back to m-space + + 1. Apply geometric matrices S,T,X,Y,Z in m-space + 2. Inverse DFT to θ-space + 3. Divide by J·B² at each θ + 4. Forward DFT back to m-space """ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_intr, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + equil::Equilibrium.PlasmaEquilibrium, + metric::ForceFreeStates.MetricData) # Copy mode numbers from FFS kf_intr.mlow = ffs_intr.mlow kf_intr.mhigh = ffs_intr.mhigh @@ -372,9 +412,9 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in psi = psi_grid[ipsi] # Get Clebsch displacement vectors at this ψ - xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] + xsp = view(xi_modes.clebsch_psi, ipsi, :) # ξ^ψ [mpert] xmp1 = view(xi_modes.clebsch_psi1, ipsi, :) # ∂ξ^ψ/∂ψ [mpert] - xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] + xms = view(clebsch_alpha_mat, ipsi, :) # ξ^α [mpert] # Evaluate geometric matrices at ψ → mpert² flat vectors, reshape to mpert×mpert geom_mats.smats(smat_flat, psi; hint=hint_s) @@ -394,8 +434,8 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state, ffs_in mul!(jbb_kapx, smat, xsp) mul!(jbb_kapx, tmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += tmat * xms mul!(jbb_divx, xmat, xmp1) - mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp - mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms + mul!(jbb_divx, ymat, xsp, 1.0 + 0.0im, 1.0 + 0.0im) # += ymat * xsp + mul!(jbb_divx, zmat, xms, 1.0 + 0.0im, 1.0 + 0.0im) # += zmat * xms @. jbb_dbob = -(jbb_divx + jbb_kapx) # Inverse DFT to θ-space, divide by J·B², forward DFT back @@ -422,9 +462,9 @@ Matches Fortran set_peq lines 859-868: transforms JBB-weighted m-space data to θ-space, removes the J·B² weighting at each poloidal angle, and transforms back. """ function _jbb_deweight!(out::AbstractVector{ComplexF64}, jbb_modes::Vector{ComplexF64}, - ft::Utilities.FourierTransforms.FourierTransform, - psi::Float64, equil::Equilibrium.PlasmaEquilibrium, - mthsurf::Int, theta_buf::Vector{ComplexF64}) + ft::Utilities.FourierTransforms.FourierTransform, + psi::Float64, equil::Equilibrium.PlasmaEquilibrium, + mthsurf::Int, theta_buf::Vector{ComplexF64}) # Inverse DFT: m-space → θ-space theta_buf .= Utilities.FourierTransforms.inverse(ft, jbb_modes) @@ -494,8 +534,8 @@ Accumulated results from all KineticForces computations. Written to gpec.h5 under the "kinetic_forces" group. """ @kwdef mutable struct KineticForcesState - method_results::Dict{String, MethodResult} = Dict{String, MethodResult}() + method_results::Dict{String,MethodResult} = Dict{String,MethodResult}() # Block-diagonal kinetic matrices: key=method, value=(numpert_total, numpert_total, 6) - kinetic_matrices::Dict{String, Array{ComplexF64,3}} = Dict{String, Array{ComplexF64,3}}() + kinetic_matrices::Dict{String,Array{ComplexF64,3}} = Dict{String,Array{ComplexF64,3}}() completed::Bool = false end diff --git a/src/KineticForces/Output.jl b/src/KineticForces/Output.jl index d1fef4e0..6fadf73a 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -15,8 +15,8 @@ Write all KineticForces results to the "kinetic_forces" group in gpec.h5. - `h5file::HDF5.File`: Open HDF5 file handle - `state::KineticForcesState`: Accumulated computation results """ -function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState) - g = create_group(h5file, "kinetic_forces") +function write_to_hdf5!(h5file::HDF5.File, state::KineticForcesState; group_name::AbstractString="kinetic_forces") + g = create_group(h5file, group_name) for (method_name, result) in state.method_results mg = create_group(g, method_name) diff --git a/src/KineticForces/Torque.jl b/src/KineticForces/Torque.jl index af22d165..17c93133 100644 --- a/src/KineticForces/Torque.jl +++ b/src/KineticForces/Torque.jl @@ -153,6 +153,8 @@ function tpsi!(tpsi_var::Ref{ComplexF64}, psi::Float64, n::Int, l::Int, dT_s_dpsi = kinetic_profiles.Te_deriv(psi) nu_s = kinetic_profiles.nue_spline(psi) else + # ni_spline / nui_spline carry THIS species' resonant density and its + # full-composition collisionality (per-species views built by the loader). n_s = kinetic_profiles.ni_spline(psi) T_s = kinetic_profiles.Ti_spline(psi) dn_s_dpsi = kinetic_profiles.ni_deriv(psi) diff --git a/test/runtests.jl b/test/runtests.jl index 9bd7b0ed..221ab8de 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -33,6 +33,7 @@ else include("./runtests_innerlayer.jl") include("./runtests_tj_analytic.jl") include("./runtests_kinetic.jl") + include("./runtests_multiion.jl") include("./runtests_fullruns.jl") include("./runtests_coils.jl") include("./runtests_imas.jl") diff --git a/test/runtests_multiion.jl b/test/runtests_multiion.jl new file mode 100644 index 00000000..d02bb7bf --- /dev/null +++ b/test/runtests_multiion.jl @@ -0,0 +1,125 @@ +using HDF5 + +@testset "Multi-ion NTV" begin + EQ = GeneralizedPerturbedEquilibrium.Equilibrium + KF = GeneralizedPerturbedEquilibrium.KineticForces + + @testset "multi_ion_composition" begin + ne = 1.24e20 + ni = 1.13e20 + zimp = 6 + mimp = 12 + # Reduces exactly to the single-ion formula Zeff = zimp - (ni/ne)*zi*(zimp-zi). + for zi in (1, 2) + zeff_old = zimp - (ni / ne) * zi * (zimp - zi) + zeff, zpitch, nmain, nimp = EQ.multi_ion_composition([zi], [ni], ne; zimp=zimp, mimp=mimp) + @test zeff ≈ zeff_old rtol = 1e-12 + @test nmain ≈ ni + end + # 50/50 D-T with total ni == single ion at total ni; quasineutrality closes the impurity. + zdt, _, nmdt, nimpdt = EQ.multi_ion_composition([1, 1], [ni / 2, ni / 2], ne; zimp=zimp, mimp=mimp) + zsingle, _, _, _ = EQ.multi_ion_composition([1], [ni], ne; zimp=zimp, mimp=mimp) + @test zdt ≈ zsingle rtol = 1e-12 + @test nmdt ≈ ni + @test ne ≈ ni + zimp * nimpdt rtol = 1e-12 # ne = Σ z_s n_s + z_imp n_imp + end + + @testset "IonSpecies + TOML conversion" begin + # TOML.jl parses [[KineticForces.ion_species]] as a Vector{Any} of Dicts. + v = Any[Dict("z" => 1, "m" => 2, "fraction" => 0.5), Dict("z" => 1, "m" => 3, "fraction" => 0.5)] + c = KF.KineticForcesControl(; ion_species=v) + @test length(c.ion_species) == 2 + @test (c.ion_species[1].z, c.ion_species[1].m, c.ion_species[1].fraction) == (1, 2, 0.5) + @test isempty(KF.KineticForcesControl().ion_species) # default single-ion + end + + @testset "resolve_ntv_species" begin + # Synthetic ASCII kinetic profile: psi, n_i(total), n_e, T_i, T_e, omega_E. + psi = collect(0.0:0.05:1.0) + Ni = @. 1.0e20 * (1 - 0.5psi) + ne = @. 1.1e20 * (1 - 0.5psi) + Ti = @. 2000.0 * (1 - 0.5psi) + Te = @. 2500.0 * (1 - 0.5psi) + wE = @. 1.0e4 * (1 - psi) + f = tempname() * ".txt" + open(f, "w") do io + println(io, "# psi ni ne Ti Te wexb") + for i in eachindex(psi) + println(io, join((psi[i], Ni[i], ne[i], Ti[i], Te[i], wE[i]), " ")) + end + end + + # Single species [z=1,m=2,fraction=1] reproduces load_kinetic_profiles. + single = EQ.load_kinetic_profiles(f; zi=1, zimp=6, mi=2, mimp=12) + sp1 = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=1.0)]; electron=false, zimp=6, mimp=12) + main = sp1[1] + for ψ in (0.3, 0.6) + @test main.profiles.ni_spline(ψ) ≈ single.ni_spline(ψ) rtol = 1e-8 + @test main.profiles.nui_spline(ψ) ≈ single.nui_spline(ψ) rtol = 1e-8 + @test main.profiles.zeff_spline(ψ) ≈ single.zeff_spline(ψ) rtol = 1e-8 + end + + # D-T-e set = {D, T, impurity, electron}; ν scales as 1/√m at fixed z, T. + sp = EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2, fraction=0.5), KF.IonSpecies(; z=1, m=3, fraction=0.5)]; + electron=true, zimp=6, mimp=12) + @test length(sp) == 4 + @test count(s -> !s.electron && s.z == 1, sp) == 2 # D, T + @test any(s -> s.z == 6 && !s.electron, sp) # impurity + @test any(s -> s.electron, sp) # electron + iD = findfirst(s -> s.m == 2, sp) + iT = findfirst(s -> s.m == 3, sp) + @test sp[iD].profiles.nui_spline(0.5) / sp[iT].profiles.nui_spline(0.5) ≈ sqrt(3 / 2) rtol = 1e-6 + + # Exactly one of {fraction, density} required. + @test_throws ErrorException EQ.resolve_ntv_species(f, [KF.IonSpecies(; z=1, m=2)]; zimp=6, mimp=12) + rm(f) + end + + @testset "explicit per-species profiles (HDF5)" begin + psi = collect(0.0:0.05:1.0) + Ni = @. 1.0e20 * (1 - 0.5psi) + ne = @. 1.1e20 * (1 - 0.5psi) + Ti = @. 2000.0 * (1 - 0.5psi) + Te = @. 2500.0 * (1 - 0.5psi) + wE = @. 1.0e4 * (1 - psi) + h5 = tempname() * ".h5" + HDF5.h5open(h5, "w") do file + file["psi"] = psi + file["n_e"] = ne + file["T_i"] = Ti + file["T_e"] = Te + file["omega_E"] = wE + file["n_i"] = Ni + file["n_D"] = 0.6 .* Ni + file["n_T"] = 0.4 .* Ni # unequal explicit shapes + end + sp = EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_D"), KF.IonSpecies(; z=1, m=3, density="n_T")]; + electron=false, zimp=6, mimp=12) + iD = findfirst(s -> s.m == 2, sp) + iT = findfirst(s -> s.m == 3, sp) + @test sp[iD].profiles.ni_spline(0.5) ≈ 0.6 * 1.0e20 * (1 - 0.25) rtol = 1e-6 + @test sp[iT].profiles.ni_spline(0.5) ≈ 0.4 * 1.0e20 * (1 - 0.25) rtol = 1e-6 + # A missing named profile errors clearly. + @test_throws ErrorException EQ.resolve_ntv_species(h5, [KF.IonSpecies(; z=1, m=2, density="n_He")]; zimp=6, mimp=12) + rm(h5) + end + + @testset "combine_species_states" begin + # Combined profile = Σ species dT/dψ interpolated onto the union grid; total = Σ total. + mk(psi, dt) = ( + st = KF.KineticForcesState(); + st.method_results["fgar"] = KF.MethodResult(; method="fgar", nn=1, + total_torque=ComplexF64(sum(dt)), psi_grid=collect(Float64, psi), + dtdpsi=ComplexF64.(dt), t_cumulative=zeros(ComplexF64, length(dt))); + st + ) + s1 = mk(0.0:0.1:1.0, fill(1.0, 11)) + s2 = mk(0.05:0.1:0.95, fill(2.0, 10)) # different (offset) grid + c = KF.combine_species_states([s1, s2]) + r = c.method_results["fgar"] + @test real(r.total_torque) ≈ 31.0 # 11·1 + 10·2 + @test maximum(abs.(r.dtdpsi)) > 0 # not the all-zeros regression + overlap = [real(d) for (p, d) in zip(r.psi_grid, r.dtdpsi) if 0.1 <= p <= 0.9] + @test all(x -> isapprox(x, 3.0; atol=1e-9), overlap) # 1 + 2 on the overlap + end +end