diff --git a/.gitignore b/.gitignore index fcb3cc96c..1d1bf86a2 100644 --- a/.gitignore +++ b/.gitignore @@ -38,3 +38,6 @@ scratch/ # Local profiling scratch (one-off study scripts, not part of the package) profiling/ + +# Override-case rundirs materialized beside their example by the regression harness +.regress-override-* diff --git a/docs/development/hdf5-conventions.md b/docs/development/hdf5-conventions.md index c6f568931..550702798 100644 --- a/docs/development/hdf5-conventions.md +++ b/docs/development/hdf5-conventions.md @@ -25,7 +25,7 @@ These rules govern `gpec.h5` (and any future GPEC-produced HDF5 output); harness - **"rational" over "singular"** in dataset names (`rational_psi`, `rational_q`, `rational_m`, `rational_n`, `rational_index`, `rational_count`) — kinetic/resistive runs are not singular at the rationals. Specifier order is standardized specifier-first (`rational_psi`, never `psi_rational`). - **Vector components** follow `[d][_]_[dpsi]` — the variable always comes first and the coordinate is always the trailing subscript. A bare coordinate suffix is the **contravariant** component (`xi_psi` = ξ^ψ), `_cov_` marks the **covariant** one (`b_cov_theta` = b_θ), a leading `J` marks a **Jacobian-weighted** component (`Jxi_theta` = J·ξ^θ), and other representations sit in the same slot (`xi_clebsch_psi`, `dxi_clebsch_psidpsi`). Never drop the variable: `clebsch_psi` is wrong because it does not say *what* is being represented. There is no HDF5/netCDF standard for super- vs subscripts — flat `_` names are universal — so the typeset form always appears in the dataset's `long_name`. - **Coordinates**: the radial abscissa is `psi` (normalized poloidal flux ψ_N) and the poloidal one is `theta` in every group; never `psi_n`, `xs`, or `ys`. -- **One name per physical quantity**: a quantity written in several groups carries the identical leaf name everywhere (`rational_psi` in `SingularSurfaces/`, `Solutions/GalerkinIntegration/`, and `SingularCoupling/`; `Delta_prime_matrix` in `SingularSurfaces/` and `Tearing/PerSurface/`; `dVdpsi` in `Profiles/`, `SingularSurfaces/`, and `KineticForces//`) — the group supplies the context, the leaf supplies the identity. +- **One name per physical quantity**: a quantity written in several groups carries the identical leaf name everywhere (`rational_psi` in `SingularSurfaces/`, `Solutions/GalerkinIntegration/`, and `SingularCoupling/`; `Delta_prime_matrix` in `SingularSurfaces/` and, on the GGJ path, `Tearing/PerSurface/` — the r_s-referenced matrix the slab layer matches is a different quantity and so carries its own name, `Tearing/PerSurface/Delta_prime_matrix_rs`; `dVdpsi` in `Profiles/`, `SingularSurfaces/`, and `KineticForces//`) — the group supplies the context, the leaf supplies the identity. ## Inputs live only under `Input/` diff --git a/docs/src/kinetic_forces.md b/docs/src/kinetic_forces.md index 858bb53a0..95cf10315 100644 --- a/docs/src/kinetic_forces.md +++ b/docs/src/kinetic_forces.md @@ -138,6 +138,35 @@ terms respectively, and do not modify the stored kinetic profile splines. when computing the `toroidal_rotation_factor` back-solve. Julia uses a clean reimplementation with consistent pre-scaling derivatives throughout. +## Regularization: why kinetic runs set `reg_spot = 0` + +`[PerturbedEquilibrium] reg_spot` smooths the displacements before they drive the NTV +integrand, multiplying ``\xi^{\psi\prime}`` and ``\xi^\alpha`` by +``Q^2/(Q^2 + \mathrm{reg\_spot}^2)`` with ``Q = m - nq``. It exists because **ideal** MHD is +singular at the rationals: ``\bar F_\mathrm{ideal} = Q F Q`` has ``\det \bar F = 0`` there, so +those two components diverge as ``1/Q`` and the torque integral does not converge. + +The **self-consistent kinetic** workflow has no such singularity. Park & Logan +([Phys. Plasmas 24, 032505 (2017)](https://doi.org/10.1063/1.4978562), §III D) decompose the +kinetic composite matrix as ``F_k = Q \bar F_k Q - P_l^\dagger Q - Q P_u + R_1`` where +``R_1 \neq 0`` at ``Q = 0``; with finite torque ``\det \bar F`` is complex, its zeros leave the +real ``\psi`` axis, and the singularity is removed from both the solution and the torque +integral. (Torque-free kinetic energy principles instead *shift and split* the zeros to +``\psi_r \mp r_{L,R}``, where the singularity is logarithmic and integrable — still not a case +for smoothing.) Logan, *Electromagnetic Torque in Tokamaks with Toroidal Asymmetries* (PhD thesis, +Princeton, 2015), Ch. 7 Eq. 7.46 states the same resonance shifting/splitting result in GPEC's own +notation. + +GPEC therefore **forces `reg_spot = 0` whenever `kinetic_factor > 0`**, logging the override. +Leaving it on suppresses a finite physical response and does so inconsistently — ``\xi^\psi`` +is never regularized, so damping the other two breaks their near-resonance cancellation in +``\delta B/B`` and leaves a spurious residue driving the NTV integrand. + +The check that this is right is that the NTV torque and the Euler–Lagrange solution's own +dissipation ``-2n\,\mathrm{Im}\langle \xi, u_2\rangle/4\mu_0`` — two independent calculations +of the same quantity — agree to a fraction of a percent with the regularization off, and disagree +by tens of percent with it on. + ## HDF5 outputs: complex torque convention and the EnergyIntegrals layout The method level of `KineticForces//` reports the two physical scalars a user diff --git a/examples/DIIID-like_SLAYER_example/gpec.toml b/examples/DIIID-like_SLAYER_example/gpec.toml index 487d41f9b..da07c4a9a 100644 --- a/examples/DIIID-like_SLAYER_example/gpec.toml +++ b/examples/DIIID-like_SLAYER_example/gpec.toml @@ -61,8 +61,11 @@ dmlim = 0.2 # Truncate integration at (last_rational_q + dmli # Runs in the force_termination path (PE skipped). Kinetic profiles (incl. the # χ⊥/χ_φ heat/momentum diffusivities) are read from the standardized HDF5 # kinetic file shared with the kinetic/NTV physics. Per-surface (uncoupled) -# analysis: the headline is the unstable 2/1; the validity gate drops surfaces -# with no real root (e.g. the 5/1, whose Δ' BVP yields a huge uncancellable Δ'). +# analysis: track the 2/1, 3/1, and 4/1 (the headline is the unstable 2/1). +# The 5/1 and 6/1 sit near the edge (ψ_N ≥ 0.97), far from the core where the +# slab-layer approximation holds; their Δ' are excessively large and likely +# invalid, no Re/Im contour crossing polishes to a true zero, and the validity +# gate correctly reports them as no_root rather than a spurious growth rate. enabled = true # Run the SLAYER tearing-mode analysis inner_model = "slayer_fitzpatrick" # Inner-layer Δ(Q) model: "slayer_fitzpatrick", "ggj_shooting", or "ggj_galerkin" scan_mode = "amr" # Q-plane scan strategy: "amr" (adaptive refinement) or "brute_force" diff --git a/examples/DIIID-like_ideal_example/gpec.toml b/examples/DIIID-like_ideal_example/gpec.toml index 24b2b7f24..9896f4c89 100644 --- a/examples/DIIID-like_ideal_example/gpec.toml +++ b/examples/DIIID-like_ideal_example/gpec.toml @@ -87,5 +87,6 @@ f0type = "maxwellian" # Equilibrium distribution moment = "pressure" # Pressure-moment NTV torque atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit write_outputs_to_HDF5 = true # Write outputs to the HDF5 file verbose = true # Enable verbose logging diff --git a/examples/Solovev_kinetic_NTV_example/gpec.toml b/examples/Solovev_kinetic_NTV_example/gpec.toml index 39483c808..0e45a65fd 100644 --- a/examples/Solovev_kinetic_NTV_example/gpec.toml +++ b/examples/Solovev_kinetic_NTV_example/gpec.toml @@ -87,3 +87,4 @@ f0fac = 1 # Scale toroidal field at constant pressure (β, q change; Φ, p, [KineticForces] kinetic_file = "kinetic.dat" # Kinetic profile file: psi_n, n_i, n_e, T_i, T_e, omega_E columns +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/examples/a10_kinetic_example/gpec.toml b/examples/a10_kinetic_example/gpec.toml index 47e2e603b..7eabffae0 100644 --- a/examples/a10_kinetic_example/gpec.toml +++ b/examples/a10_kinetic_example/gpec.toml @@ -59,3 +59,4 @@ nutype = "harmonic" # Collision operator (zero, small, krook, harmoni f0type = "maxwellian" # Distribution function (maxwellian, jkp, cgl) atol_xlmda = 1e-9 # Absolute tolerance for inner pitch + energy integrations rtol_xlmda = 1e-5 # Relative tolerance for inner pitch + energy integrations +axis_validity_suppression = true # Suppress kinetic terms inside the near-axis orbit-width validity limit diff --git a/regression-harness/cases/diiid_n1.toml b/regression-harness/cases/diiid_n1.toml index feebc0c7f..dcba09788 100644 --- a/regression-harness/cases/diiid_n1.toml +++ b/regression-harness/cases/diiid_n1.toml @@ -411,3 +411,14 @@ extract = "value" label = "Runtime (s)" noise_threshold = 0.0 order = 999 + +# Root-area-weighted (coordinate-invariant) forcing field b̃ on the control surface: the +# applied coil spectrum pushed through the √(J|∇ψ|) weight operator. Tracks the field-space +# conversion that every ResponseMatrices/ field-space output depends on. +[quantities.forcing_b_rootarea_norm] +h5path = "PerturbedEquilibrium/forcing_b_root_area" +type = "complex_vector" +extract = "norm" +label = "||forcing b~|| (root-area-weighted)" +noise_threshold = 1e-10 +order = 1000 diff --git a/regression-harness/cases/diiid_slayer_n1.toml b/regression-harness/cases/diiid_slayer_n1.toml index 540ca7e50..90edf2e09 100644 --- a/regression-harness/cases/diiid_slayer_n1.toml +++ b/regression-harness/cases/diiid_slayer_n1.toml @@ -1,6 +1,8 @@ # Regression case: DIII-D-like H-mode n=1 SLAYER tearing-mode growth rates. # Uncoupled per-surface analysis with the AMR Q-plane scan; tracks the per-surface -# layer parameters and the extracted roots (Q, ω, γ) for the inner three rationals. +# layer parameters, the diamagnetic drive Q_e/Q_i, and the extracted roots (Q, ω, γ) +# for the inner three rationals. The n-scaling of Q_e/Q_i, which no single-n case can +# constrain, is covered by a unit test in test/runtests_slayer_inputs.jl instead. # 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] @@ -89,6 +91,24 @@ label = "SLAYER iota_e" noise_threshold = 1e-12 order = 19 +# Normalized diamagnetic frequencies. iota_e above is their ratio, so it is blind +# to a common factor on both; these pin the drive itself. +[quantities.slayer_Q_e] +h5path = "Tearing/PerSurface/Q_e" +type = "real_vector" +extract = "all_real" +label = "SLAYER Q_e" +noise_threshold = 1e-12 +order = 20 + +[quantities.slayer_Q_i] +h5path = "Tearing/PerSurface/Q_i" +type = "real_vector" +extract = "all_real" +label = "SLAYER Q_i" +noise_threshold = 1e-12 +order = 21 + # Tearing eigenvalue (coupled mode → length 1). The headline deliverable: # a real, nonzero growth rate on a realistic equilibrium. Root-extraction # is sensitive to the AMR cell topology and ODE solver, so the Q/ω/γ diff --git a/regression-harness/src/runner.jl b/regression-harness/src/runner.jl index ed7e5f4d5..1a192a140 100644 --- a/regression-harness/src/runner.jl +++ b/regression-harness/src/runner.jl @@ -27,16 +27,18 @@ end Materialize the directory GPEC will actually run in. With no `overrides`, the example deck runs in place (returns it untouched). With overrides, -the deck is copied to a throwaway temp dir and the named gpec.toml keys are patched there, -so one shared example can serve several cases (e.g. a collisionless variant via -`"KineticForces.nutype" => "zero"`). Override keys are dotted paths into the TOML; missing -intermediate tables are created. Returns `(rundir, is_temp)`; the caller removes the temp -tree when `is_temp`. +the deck is copied to a uniquely named throwaway sibling of the example dir and the named gpec.toml keys are +patched there, so one shared example can serve several cases (e.g. a collisionless variant via +`"KineticForces.nutype" => "zero"`). The copy sits beside the original so relative file +references in the deck (e.g. `SLAYER.profile_file = "..//..."`) still resolve. +Override keys are dotted paths into the TOML; missing intermediate tables are created. +Returns `(rundir, is_temp)`; the caller removes `rundir` when `is_temp`. """ function _materialize_rundir(example_path::String, overrides::Dict{String,Any}) isempty(overrides) && return (example_path, false) - rundir = joinpath(mktempdir(), "deck") - cp(example_path, rundir) + # Unique per invocation, so concurrent runs from one checkout cannot share (and race on) a rundir. + rundir = mktempdir(dirname(example_path); prefix=".regress-override-" * basename(example_path) * "-", cleanup=false) + cp(example_path, rundir; force=true) rm(joinpath(rundir, "gpec.h5"); force=true) # drop any stale output copied along cfg = TOML.parsefile(joinpath(rundir, "gpec.toml")) for (dotted, val) in overrides @@ -523,7 +525,8 @@ function run_local(db::SQLite.DB, case_spec::CaseSpec, repo_root::String; rm(runinfo_file; force=true) end if rundir_is_temp && rundir !== nothing - rm(dirname(rundir); recursive=true, force=true) + # rundir sits beside the example; removing its parent would delete examples/ itself + rm(rundir; recursive=true, force=true) end end end @@ -647,7 +650,8 @@ function run_at_commit(db::SQLite.DB, commit_hash::String, ref_name::String, rm(runinfo_file; force=true) end if rundir_is_temp && rundir !== nothing - rm(dirname(rundir); recursive=true, force=true) + # rundir sits beside the example; removing its parent would delete examples/ itself + rm(rundir; recursive=true, force=true) end if own_worktree && worktree_path !== nothing remove_worktree(worktree_path, repo_root) diff --git a/src/Analysis/CoilForcing.jl b/src/Analysis/CoilForcing.jl index 050b8cf45..660d0505b 100644 --- a/src/Analysis/CoilForcing.jl +++ b/src/Analysis/CoilForcing.jl @@ -238,11 +238,7 @@ function plot_mode_spectrum(forcing_modes; mlow=nothing, mhigh=nothing, save_pat end # Helicity = sign(Bt)·sign(Ip); the SFL→machine toroidal-angle handedness. -function _equil_helicity(equil) - bt = isnothing(equil.params.bt_sign) ? 1 : Int(sign(equil.params.bt_sign)) - ip = isnothing(equil.params.crnt) ? 1 : Int(sign(equil.params.crnt)) - return bt * ip -end +_equil_helicity(equil) = equil.params.bt_sign * equil.params.ip_sign # Unpack `modes` (a Vector of ForcingMode-like objects filtered to toroidal mode `n`, # or a `(m_vals, amplitudes)` tuple) into parallel m and complex-amplitude vectors. diff --git a/src/Equilibrium/AnalyticEquilibrium.jl b/src/Equilibrium/AnalyticEquilibrium.jl index 72a16a05c..767f2b718 100644 --- a/src/Equilibrium/AnalyticEquilibrium.jl +++ b/src/Equilibrium/AnalyticEquilibrium.jl @@ -242,10 +242,10 @@ Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ function tj_analytic_f1(x::Float64, nu::Float64, qc::Float64) if x < 0.1 x2 = x * x - return x2 * (1 - (nu-1)*x2/2 + (nu-1)*(nu-2)*x2*x2/6 - - (nu-1)*(nu-2)*(nu-3)*x2*x2*x2/24) / qc + return x2 * (1 - (nu - 1) * x2 / 2 + (nu - 1) * (nu - 2) * x2 * x2 / 6 - + (nu - 1) * (nu - 2) * (nu - 3) * x2 * x2 * x2 / 24) / qc else - return (1 - (1 - x*x)^nu) / (nu * qc) + return (1 - (1 - x * x)^nu) / (nu * qc) end end @@ -259,10 +259,10 @@ parameterization. function tj_analytic_f1p(x::Float64, nu::Float64, qc::Float64) if x < 0.1 x2 = x * x - return 2*x * (1 - (nu-1)*x2 + (nu-1)*(nu-2)*x2*x2/2 - - (nu-1)*(nu-2)*(nu-3)*x2*x2*x2/6) / qc + return 2 * x * (1 - (nu - 1) * x2 + (nu - 1) * (nu - 2) * x2 * x2 / 2 - + (nu - 1) * (nu - 2) * (nu - 3) * x2 * x2 * x2 / 6) / qc else - return 2*x * (1 - x*x)^(nu-1) / qc + return 2 * x * (1 - x * x)^(nu - 1) / qc end end @@ -273,10 +273,11 @@ GPEC adaptation of the analytic shape ODE used in R. Fitzpatrick's TJ code `tj_analytic_run_direct` call so both pipelines share identical numerics. Fields: + - physical: a, R0, qc, mu, pc, B0 - derived: epsa2 = (a/R0)² - near-axis BC constants: rmin, x0 = rmin, r0 = rmin·a, f1c = 1/qc, - p2ppc = d²p₂/dx²|_0 = −2·μ·pc + p2ppc = d²p₂/dx²|_0 = −2·μ·pc """ struct TJAnalyticShapeParams a::Float64 @@ -293,15 +294,15 @@ struct TJAnalyticShapeParams p2ppc::Float64 end -function TJAnalyticShapeParams(tj::TJAnalyticConfig; rmin::Float64 = 1e-4) +function TJAnalyticShapeParams(tj::TJAnalyticConfig; rmin::Float64=1e-4) a, R0 = tj.lar_a, tj.lar_r0 - mu = max(tj.mu, 1.001) + mu = max(tj.mu, 1.001) return TJAnalyticShapeParams( a, R0, tj.qc, mu, tj.pc, tj.B0, (a / R0)^2, rmin, rmin, rmin * a, 1.0 / tj.qc, - -2.0 * mu * tj.pc, + -2.0 * mu * tj.pc ) end @@ -315,9 +316,9 @@ The params argument carries TJAnalyticShapeParams fields plus the current `nu`. """ function tj_analytic_shape_rhs!(dy, y, params, r) (; a, B0, qc, mu, pc, epsa2, nu) = params - x = r / a + x = r / a xfac = max(1 - x^2, 0.0) - f1 = tj_analytic_f1(x, nu, qc) + f1 = tj_analytic_f1(x, nu, qc) f1px = tj_analytic_f1p(x, nu, qc) p2px = -2 * mu * pc * x * xfac^(mu - 1) @@ -340,16 +341,19 @@ function tj_analytic_shape_rhs!(dy, y, params, r) # f₃'(x) for Hₙ = Vₙ = 0 (n ≥ 2 harmonics rescaled to zero, as in the # TJ-analytic benchmark configuration of Fitzpatrick's TJ code). g2, f3 = y[2], y[5] - f3p_x = -f3 * f1px / f1 - - f1 * (3 * x^2 / 2 - 2 * x * H1p + H1p^2) / x + - f1px * (g2 - 3 * x^2 / 4 + H1 + 3 * H1p^2 / 2) + - x^2 * p2px * (g2 + x^2 / 2 - 3 * x * H1p - 2 * H1) / f1 + f3p_x = + -f3 * f1px / f1 - + f1 * (3 * x^2 / 2 - 2 * x * H1p + H1p^2) / x + + f1px * (g2 - 3 * x^2 / 4 + H1 + 3 * H1p^2 / 2) + + x^2 * p2px * (g2 + x^2 / 2 - 3 * x * H1p - 2 * H1) / f1 dy[5] = f3p_x / a return nothing end -"""Initial conditions at x = x0, matching the TJ-analytic model's near-axis -expansion (cf. R. Fitzpatrick's TJ code, https://github.com/rfitzp/TJ).""" +""" +Initial conditions at x = x0, matching the TJ-analytic model's near-axis +expansion (cf. R. Fitzpatrick's TJ code, https://github.com/rfitzp/TJ). +""" function tj_analytic_shape_initial(p::TJAnalyticShapeParams, nu::Float64) f1_0 = tj_analytic_f1(p.x0, nu, p.qc) y0 = zeros(5) @@ -368,14 +372,14 @@ downstream Hₙ / ψ splines sit on uniform nodes); leave it `nothing` for the default adaptive save pattern used by `tj_analytic_run`. """ function tj_analytic_shape_solve(p::TJAnalyticShapeParams, nu::Float64; - reltol::Float64 = 1e-7, abstol::Float64 = 1e-8, - saveat = nothing) - rhs_params = (; p.a, p.B0, p.qc, p.mu, p.pc, p.epsa2, nu = nu) + reltol::Float64=1e-7, abstol::Float64=1e-8, + saveat=nothing) + rhs_params = (; p.a, p.B0, p.qc, p.mu, p.pc, p.epsa2, nu=nu) prob = ODEProblem(tj_analytic_shape_rhs!, tj_analytic_shape_initial(p, nu), (p.r0, p.a), rhs_params) if saveat === nothing - return solve(prob, Vern9(); reltol, abstol, maxiters = 10000, dense = false) + return solve(prob, Vern9(); reltol, abstol, maxiters=10000, dense=false) else - return solve(prob, Vern9(); reltol, abstol, maxiters = 10000, saveat = saveat) + return solve(prob, Vern9(); reltol, abstol, maxiters=10000, saveat=saveat) end end @@ -389,9 +393,9 @@ O(εa²) correction relative to the lowest-order guess ν = qa/qc, which matters for the TJ-analytic benchmark at large ε. Falls back to the lowest-order ν if the bracket search diverges. """ -function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; reltol::Float64 = 1e-7) +function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; reltol::Float64=1e-7) function q2_edge(nu::Float64) - sol = tj_analytic_shape_solve(p, nu; reltol) + sol = tj_analytic_shape_solve(p, nu; reltol) g2end = sol.u[end][2] f3end = sol.u[end][5] f1end = tj_analytic_f1(1.0, nu, p.qc) @@ -400,7 +404,7 @@ function tj_analytic_find_nu(p::TJAnalyticShapeParams, qa_target::Float64; relto nu_guess = qa_target / p.qc return try find_zero(nu -> q2_edge(nu) - qa_target, (0.5 * nu_guess, 2 * nu_guess); - atol = 1e-8, rtol = 1e-10) + atol=1e-8, rtol=1e-10) catch err @warn "ν root-find failed for TJ-analytic equilibrium; falling back to lowest-order ν = qa/qc" error = err nu_guess @@ -443,16 +447,16 @@ included; they are zero in the TJ-analytic benchmark scans. Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ """ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig) - a, R0 = tj.lar_a, tj.lar_r0 + a, R0 = tj.lar_a, tj.lar_r0 qc, mu = tj.qc, max(tj.mu, 1.001) pc, B0 = tj.pc, tj.B0 ma, mtau = tj.ma, tj.mtau p = TJAnalyticShapeParams(tj) - epsa2 = p.epsa2 - p00_phys = B0^2 * epsa2 * pc # μ₀P = B₀²·εa²·p₂ at axis + epsa2 = p.epsa2 + p00_phys = B0^2 * epsa2 * pc # μ₀P = B₀²·εa²·p₂ at axis - nu = tj_analytic_find_nu(p, tj.qa; reltol = equil_input.etol) - sol = tj_analytic_shape_solve(p, nu; reltol = equil_input.etol) + nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol) + sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol) r_arr = sol.t y_mat = reduce(hcat, sol.u)' @@ -468,7 +472,7 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig) xfac = max(1 - x^2, 0.0) f1 = tj_analytic_f1(x, nu, qc) - ψ = y_mat[i, 1] + ψ = y_mat[i, 1] g2 = y_mat[i, 2] H1 = y_mat[i, 3] f3 = y_mat[i, 5] @@ -507,7 +511,7 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig) f = spl(r; hint=hint) # f[1]=F, f[2]=P, f[3]=q, f[4]=ψ, f[5]=g₂, f[6]=H₁ - sq_xs[ia] = f[4] / psio + sq_xs[ia] = f[4] / psio sq_fs[ia, 1] = f[1] # F sq_fs[ia, 2] = f[2] # P sq_fs[ia, 3] = f[3] # q @@ -525,7 +529,7 @@ function tj_analytic_run(equil_input::EquilibriumConfig, tj::TJAnalyticConfig) for itau in 1:(mtau+1) θ = 2π * (itau - 1) / mtau rzphi_fs_nodes[ia, itau, 1] = R0 + Δ + α * r * cos(θ) - rzphi_fs_nodes[ia, itau, 2] = α * r * sin(θ) + rzphi_fs_nodes[ia, itau, 2] = α * r * sin(θ) end end @@ -579,38 +583,38 @@ ODE (g₂, H₁, H₁', f₃), the `GetPSIvac` / `GetHHvac` vacuum extension, an EFIT-writer (R, Z) → (r, w) Newton inversion that this routine adapts. """ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticConfig; - nrbox::Int = 257, nzbox::Int = 257, rc::Float64 = 1.2) - a, R0 = tj.lar_a, tj.lar_r0 + nrbox::Int=257, nzbox::Int=257, rc::Float64=1.2) + a, R0 = tj.lar_a, tj.lar_r0 qc, mu = tj.qc, max(tj.mu, 1.001) pc, B0 = tj.pc, tj.B0 p = TJAnalyticShapeParams(tj) epsa, epsa2 = p.a / p.R0, p.epsa2 - p00_phys = B0^2 * epsa2 * pc + p00_phys = B0^2 * epsa2 * pc # ν root-find (cf. Fitzpatrick TJ's Setnu): q₂(1) = qa_target. - nu = tj_analytic_find_nu(p, tj.qa; reltol = equil_input.etol) + nu = tj_analytic_find_nu(p, tj.qa; reltol=equil_input.etol) # Dense saveat so the downstream splines (H₁, g₂, f₃, ψ) are evaluated on # a fine uniform r grid rather than the ~30 adaptive Vern9 steps — otherwise # the (R, Z) → (r, w) Newton iteration hits spline interpolation artifacts. - dense_r = collect(range(p.r0, p.a; length = 1024)) - sol = tj_analytic_shape_solve(p, nu; reltol = equil_input.etol, - abstol = 1e-10, saveat = dense_r) - r_arr = sol.t - y_mat = reduce(hcat, sol.u)' + dense_r = collect(range(p.r0, p.a; length=1024)) + sol = tj_analytic_shape_solve(p, nu; reltol=equil_input.etol, + abstol=1e-10, saveat=dense_r) + r_arr = sol.t + y_mat = reduce(hcat, sol.u)' # Radial splines in the TJ-analytic dimensionless x = r/a on a clean grid for H₁ etc. x_nodes = r_arr ./ a - ψ_of_r = cubic_interp(r_arr, y_mat[:, 1]; extrap=ExtendExtrap()) - H1_of_x = cubic_interp(x_nodes, y_mat[:, 3]; extrap=ExtendExtrap()) + ψ_of_r = cubic_interp(r_arr, y_mat[:, 1]; extrap=ExtendExtrap()) + H1_of_x = cubic_interp(x_nodes, y_mat[:, 3]; extrap=ExtendExtrap()) H1p_of_x = cubic_interp(x_nodes, y_mat[:, 4]; extrap=ExtendExtrap()) - g2_of_x = cubic_interp(x_nodes, y_mat[:, 2]; extrap=ExtendExtrap()) - f3_of_x = cubic_interp(x_nodes, y_mat[:, 5]; extrap=ExtendExtrap()) + g2_of_x = cubic_interp(x_nodes, y_mat[:, 2]; extrap=ExtendExtrap()) + f3_of_x = cubic_interp(x_nodes, y_mat[:, 5]; extrap=ExtendExtrap()) # Edge values needed by GetPSIvac - f1a = tj_analytic_f1(1.0, nu, qc) - f3a = f3_of_x(1.0) - H1a = H1_of_x(1.0) + f1a = tj_analytic_f1(1.0, nu, qc) + f3a = f3_of_x(1.0) + H1a = H1_of_x(1.0) H1ap = H1p_of_x(1.0) psio = ψ_of_r(a) # ψ at r = a (boundary) @@ -644,7 +648,7 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo return f_R_shift(rc - 1e-8, w) * r^2 / rc^2 end H1 = (r < 1.0) ? H1_of_x(r) : H1_vac(r) - L = r^3 / 8 - r * H1 / 2 + L = r^3 / 8 - r * H1 / 2 return epsa2 * H1 + epsa2 * epsa * L * cos(w) end function f_Z_shift(r::Float64, w::Float64) @@ -652,7 +656,7 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo return f_Z_shift(rc - 1e-8, w) * r^2 / rc^2 end H1 = (r < 1.0) ? H1_of_x(r) : H1_vac(r) - L = r^3 / 8 - r * H1 / 2 + L = r^3 / 8 - r * H1 / 2 return -epsa2 * epsa * L * sin(w) end @@ -699,9 +703,9 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo # Grid spans R₀ ± rc·a × ±rc·a (where rc is the vacuum-shell radius in # units of a), giving a comfortable margin for the separatrix finder. r_span = rc * a - psi_in_xs = collect(range(R0 - r_span, R0 + r_span; length = nrbox)) - psi_in_ys = collect(range(-r_span, r_span; length = nzbox)) - psi_rz = zeros(Float64, nrbox, nzbox) + psi_in_xs = collect(range(R0 - r_span, R0 + r_span; length=nrbox)) + psi_in_ys = collect(range(-r_span, r_span; length=nzbox)) + psi_rz = zeros(Float64, nrbox, nzbox) for i in 1:nrbox, j in 1:nzbox R_norm = psi_in_xs[i] / R0 @@ -724,8 +728,10 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo # 1D profile spline, same layout as read_efit (4 columns). Use the # TJ-analytic q₂ on the radial grid so that the prescribed q is # consistent with the ψ(R,Z) we just constructed. - psi_norm_grid = range(0.0, 1.0; length = nrbox) - F_nodes = zeros(nrbox); P_nodes = zeros(nrbox); q_nodes = zeros(nrbox) + psi_norm_grid = range(0.0, 1.0; length=nrbox) + F_nodes = zeros(nrbox) + P_nodes = zeros(nrbox) + q_nodes = zeros(nrbox) for i in 1:nrbox ψN = psi_norm_grid[i] # Invert ψN = (ψ_plasma(r) - 0) / psio ⇒ find r such that ψ_plasma(r) = ψN·psio. @@ -746,7 +752,7 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo F_nodes[i] = R0 * B0 * (1 + epsa2 * g2_val) P_nodes[i] = p00_phys * xfac^mu q_nodes[i] = (x > 1e-10) ? x^2 * (1 + epsa2 * g2_val) * - exp(-epsa2 * f3_val / f1) / f1 : qc + exp(-epsa2 * f3_val / f1) / f1 : qc end sq_fs_nodes = hcat(F_nodes, P_nodes, q_nodes, sqrt.(collect(psi_norm_grid))) sq_in = cubic_interp(collect(psi_norm_grid), Series(sq_fs_nodes); extrap=ExtendExtrap()) @@ -756,7 +762,7 @@ function tj_analytic_run_direct(equil_input::EquilibriumConfig, tj::TJAnalyticCo # Analytic: ingest=nothing — replay regenerates from the [TJ_ANALYTIC_INPUT] TOML section. return DirectRunInput(equil_input, sq_in, psi_in, psi_in_xs, psi_in_ys, - rmin_grid, rmax_grid, zmin_grid, zmax_grid, psio, 1, nothing) + rmin_grid, rmax_grid, zmin_grid, zmax_grid, psio, 1, 1, nothing) end """ @@ -836,7 +842,6 @@ function sol_run(equil_inputs::EquilibriumConfig, sol_inputs::SolovevConfig) @info "Generating Solovev equilibrium: mr=$mr, mz=$mz, ma=$ma, e=$(@sprintf("%.3f", e)), a=$(@sprintf("%.3f", a)), r0=$(@sprintf("%.3f", r0)), q0=$(@sprintf("%.3f", q0))" # 1 is bt_sign=+1: Solovev has positive Bt by construction (f0 = r0 * b0fac > 0). - # No ip_sign field is needed; sign(crnt) is recovered from params.crnt downstream. - # Analytic: ingest=nothing — replay regenerates from the [SOL_INPUT] TOML section. - return DirectRunInput(equil_inputs, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, 1, nothing) + # Analytic equilibria carry positive field and current; ingest=nothing — replay regenerates from the [SOL_INPUT] TOML section. + return DirectRunInput(equil_inputs, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, 1, 1, nothing) end diff --git a/src/Equilibrium/CoordinateInvariant.jl b/src/Equilibrium/CoordinateInvariant.jl index 497c365fd..ae001b5d6 100644 --- a/src/Equilibrium/CoordinateInvariant.jl +++ b/src/Equilibrium/CoordinateInvariant.jl @@ -54,13 +54,14 @@ w(θ) = √(J·|∇ψ|). Operationally, sqrtamat is the mode-space √weight operator: for a field b with Fourier coefficients b_fft, it satisfies the identity -`‖sqrtamat·b_fft‖² = N² · ∫ |b|² · J|∇ψ| dθ` -which is Jacobian-invariant on a given flux surface (see -`scripts/test_power_norm_invariance.jl`). - -The backward step uses exp(−imθ)/(1/N) normalization paired with the Julia -forward FT exp(+imθ) so that round-trip = identity and the convolution -structure is correct. +`‖sqrtamat·b_fft‖² = ∫ |b|² · J|∇ψ| dθ` (θ normalized to [0, 1)) +which is Jacobian-invariant on a given flux surface. Its diagonal is the θ-average of +√(J|∇ψ|), so `sqrtamat/√jarea` has a diagonal close to (and never above) one; this is the +matrix Fortran GPEC writes as `J_surf_2` after that division. + +Each column is the unit mode e^{+i m_k θ} taken to θ-space with `FourierTransforms.inverse`, +weighted pointwise, and brought back with the forward transform `ft(...)`, the same pair whose +round trip is the identity, so the convolution structure is correct. """ function compute_sqrtamat( equil::PlasmaEquilibrium, @@ -78,8 +79,8 @@ function compute_sqrtamat( e_k .= 0.0 e_k[k] = 1.0 + 0.0im - # Standard backward FT: f(θ_j) = (1/N) Σ_m c_m exp(-imθ_j) = (transpose(basis) * c) / N - theta_vec = (transpose(ft.basis) * e_k) ./ mtheta + # Inverse FT of the unit mode, f(θ_j) = exp(+i m_k θ_j), through the library's own inverse + theta_vec = Utilities.FourierTransforms.inverse(ft, e_k) # Multiply pointwise by √(J·|∇ψ|) in theta-space theta_vec .*= sqrt_jdp diff --git a/src/Equilibrium/DirectEquilibrium.jl b/src/Equilibrium/DirectEquilibrium.jl index 586d81001..512f74c54 100644 --- a/src/Equilibrium/DirectEquilibrium.jl +++ b/src/Equilibrium/DirectEquilibrium.jl @@ -706,6 +706,7 @@ robustness. params = EquilibriumParameters() params.bt_sign = raw_profile.bt_sign + params.ip_sign = raw_profile.ip_sign params.psihigh_resolved = psihigh return PlasmaEquilibrium(raw_profile.config, params, profiles, geometry, diff --git a/src/Equilibrium/EquilibriumTypes.jl b/src/Equilibrium/EquilibriumTypes.jl index 344f5a60f..ac5a212fc 100644 --- a/src/Equilibrium/EquilibriumTypes.jl +++ b/src/Equilibrium/EquilibriumTypes.jl @@ -245,7 +245,9 @@ A mutable struct holding parameters for the Large Aspect Ratio (LAR) plasma equi zeroth::Bool = false end -"Build a `LargeAspectRatioConfig` from a parsed `[LAR_INPUT]` TOML table." +""" +Build a `LargeAspectRatioConfig` from a parsed `[LAR_INPUT]` TOML table. +""" function LargeAspectRatioConfig(input_dict::Dict{String,Any}) return LargeAspectRatioConfig(; symbolize_keys(input_dict)...) end @@ -285,7 +287,9 @@ Reference: R. Fitzpatrick, TJ code, https://github.com/rfitzp/TJ zeroth::Bool = false # If true, suppress Shafranov shift end -"Build a `TJAnalyticConfig` from a parsed `[TJ_ANALYTIC_INPUT]` TOML table." +""" +Build a `TJAnalyticConfig` from a parsed `[TJ_ANALYTIC_INPUT]` TOML table. +""" function TJAnalyticConfig(input_dict::Dict{String,Any}) return TJAnalyticConfig(; symbolize_keys(input_dict)...) end @@ -321,7 +325,9 @@ A mutable struct holding parameters for the Solev'ev (SOL) plasma equilibrium mo f0fac::Float64 = 1 # scale toroidal field at constant pressure (s*f. beta,q changes. Phi,p,bp constant) end -"Build a `SolovevConfig` from a parsed `[SOL_INPUT]` TOML table." +""" +Build a `SolovevConfig` from a parsed `[SOL_INPUT]` TOML table. +""" function SolovevConfig(input_dict::Dict{String,Any}) return SolovevConfig(; symbolize_keys(input_dict)...) end @@ -346,6 +352,7 @@ not serializable; they are reconstructed from these nodes by `build_direct_from_ - `rmin/rmax/zmin/zmax::Float64` — computational-grid bounds [m] - `psio::Float64` — total flux difference |ψ_axis - ψ_boundary| [Wb/rad] - `bt_sign::Int` — sign of the toroidal field (+1 or -1) + - `ip_sign::Int` — sign of the plasma current (+1 or -1) as the source file states it """ struct DirectIngest sq_xs::Vector{Float64} @@ -359,6 +366,7 @@ struct DirectIngest zmax::Float64 psio::Float64 bt_sign::Int + ip_sign::Int end """ @@ -438,6 +446,10 @@ raw equilibrium data and preparing the initial splines. - `zmax::Float64` — Maximum Z-coordinate of the computational grid [m] - `psio::Float64` — Total flux difference `|ψ_axis - ψ_boundary|` [Wb/rad] - `bt_sign::Int` — Sign of the toroidal field (+1 or -1); read from fpol sign in EFIT g-files + - `ip_sign::Int` — Sign of the plasma current (+1 or -1); read from the `current` header value + in EFIT g-files and `global_quantities.ip` in IMAS. The internal flux is always made + positive at the axis, so the sign is kept here (never recovered from the computed current) + and fixes the SFL→machine toroidal-angle handedness `helicity = bt_sign × ip_sign`. - `ingest::EquilibriumIngest` — captured raw arrays for the `gpec.h5` rerun snapshot (a [`DirectIngest`](@ref) for file-based reads, or `nothing` for analytic equilibria) - `psihigh_resolved::Float64` — outer flux limit the equilibrium is formed on: `config.psihigh` @@ -459,6 +471,7 @@ mutable struct DirectRunInput{S<:FastInterpolations.CubicSeriesInterpolant,I2D<: zmax::Float64 # Maximum Z-coordinate of the computational grid [m]. psio::Float64 # The total flux difference |ψ_axis - ψ_boundary| [Weber / radian]. bt_sign::Int # Sign of the toroidal field: +1 or -1 (from fpol sign in g-file) + ip_sign::Int # Sign of the plasma current: +1 or -1 (from the g-file current / IMAS ip) ingest::EquilibriumIngest psihigh_resolved::Float64 end @@ -466,9 +479,9 @@ end # Readers construct without a resolved psihigh; it starts at the request and `resolve_psihigh!` # clamps it for efit-family equilibria. DirectRunInput(config::EquilibriumConfig, sq_in, psi_in, psi_in_xs, psi_in_ys, - rmin, rmax, zmin, zmax, psio, bt_sign, ingest) = + rmin, rmax, zmin, zmax, psio, bt_sign, ip_sign, ingest) = DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, - rmin, rmax, zmin, zmax, psio, bt_sign, ingest, config.psihigh) + rmin, rmax, zmin, zmax, psio, bt_sign, ip_sign, ingest, config.psihigh) """ InverseRunInput(...) @@ -610,6 +623,7 @@ A mutable struct containing computed equilibrium parameters and diagnostic flags bt0::Union{Nothing,Float64} = nothing # Toroidal magnetic field at the axis [T] (always positive; sign in bt_sign) crnt::Union{Nothing,Float64} = nothing # Plasma current at the axis [A] bt_sign::Int = 1 # Sign of the toroidal field: +1 (positive Bt) or -1 (negative Bt, e.g. DIII-D standard) + ip_sign::Int = 1 # Sign of the plasma current as the source file states it: +1 or -1 (crnt is always positive) bwall::Union{Nothing,Float64} = nothing # Toroidal magnetic field at the wall [T] verbose::Bool = false # Whether to print verbose output diagnose_src::Bool = false # Whether to diagnose source data @@ -749,7 +763,7 @@ function GeometryProfileSplines(xs::Vector{Float64}, GeometryProfileSplines{typeof(area_spline)}( xs, npts, npts - 1, - area_spline, avg_r_spline, avg_R_spline, + area_spline, avg_r_spline, avg_R_spline ) end @@ -828,7 +842,7 @@ function KineticProfileSplines(xs::Vector{Float64}, xs, npts, npts - 1, ni_spline, ne_spline, Ti_spline, Te_spline, omegaE_spline, loglam_spline, nui_spline, nue_spline, zeff_spline, - ni_deriv, ne_deriv, Ti_deriv, Te_deriv, + ni_deriv, ne_deriv, Ti_deriv, Te_deriv ) end diff --git a/src/Equilibrium/GridRefinement.jl b/src/Equilibrium/GridRefinement.jl index e58eb81dc..0b64337a7 100644 --- a/src/Equilibrium/GridRefinement.jl +++ b/src/Equilibrium/GridRefinement.jl @@ -419,6 +419,7 @@ density (`_knot_density`), equidistribution, a global minimum-spacing floor (`enforce_min_spacing`), and rational-surface bracketing (`bracket_mandatory_nodes`). `tau` is the target interpolation accuracy (`psi_accuracy`); `kin` optionally supplies kinetic profiles whose pedestal gradients attract knots; `mandatory` lists rational-surface ψ values to bracket; +`pinned` lists ψ values inserted as plain knots without a cleared zone (kinetic-resonance surfaces); `singfac_min` and `n_min` (smallest |n| in the run) set each surface's matching half-stencil `dpsi = singfac_min/(n_min·|q′|)`, and the bracket half-width is `bracket_coef·dpsi` (floored at `min_spacing`). Rational surfaces are bracketed, not pinned: a knot on the surface would make the @@ -428,6 +429,7 @@ function refined_psi_grid(equil::PlasmaEquilibrium; tau::Float64, kin::Union{Nothing,KineticProfileSplines}=nothing, mandatory::Vector{Float64}=Float64[], + pinned::Vector{Float64}=Float64[], singfac_min::Float64=1e-4, n_min::Int=1, bracket_coef::Float64=BRACKET_COEF, @@ -449,6 +451,9 @@ function refined_psi_grid(equil::PlasmaEquilibrium; N == N_cap && M_total > N_cap && @warn "refined_psi_grid: knot count capped at $N_cap (density integral wants $(ceil(Int, M_total))); psi_accuracy=$tau may not be attainable" grid = enforce_min_spacing(_equidistribute(xs, rho, N), min_spacing) + # Pinned knots use knot-at-node semantics (no cleared zone); inserted before rational bracketing + # so a rational's bracket still wins locally where the two overlap. + grid = isempty(pinned) ? grid : merge_mandatory_nodes(grid, pinned) isempty(mandatory) && return grid min_half_widths = [max(bracket_coef * singfac_min / (n_min * abs(equil.profiles.q_deriv(m))), min_spacing) for m in mandatory] return bracket_mandatory_nodes(grid, mandatory, min_half_widths, min_spacing) diff --git a/src/Equilibrium/ReadEquilibrium.jl b/src/Equilibrium/ReadEquilibrium.jl index 8cc221b7a..755a91aca 100644 --- a/src/Equilibrium/ReadEquilibrium.jl +++ b/src/Equilibrium/ReadEquilibrium.jl @@ -67,6 +67,8 @@ function read_efit(config::EquilibriumConfig) header_vals = _read_1d_gfile_format(lines[2:5], 20) rdim, zdim, rcentr, rleft, zmid = header_vals[1:5] rmaxis, zmaxis, simag, sibry = header_vals[6:9] + ip_sign = Int(sign(header_vals[11])) # the g-file's plasma current; its sign sets the helicity + ip_sign == 0 && (ip_sign = 1) # --- Parse Data Blocks --- current_line_idx = 6 @@ -120,9 +122,9 @@ function read_efit(config::EquilibriumConfig) # Capture the raw arrays that reconstruct sq_in and psi_in, so the rerun path # (gpec.h5 → setup_equilibrium) can skip the g-file parse entirely. - ingest = DirectIngest(sq_xs, sq_fs_nodes, psi_in_xs, psi_in_ys, psi_proc, rmin, rmax, zmin, zmax, psio, fpol_sign) + ingest = DirectIngest(sq_xs, sq_fs_nodes, psi_in_xs, psi_in_ys, psi_proc, rmin, rmax, zmin, zmax, psio, fpol_sign, ip_sign) - return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, fpol_sign, ingest) + return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, fpol_sign, ip_sign, ingest) end @@ -152,7 +154,7 @@ function read_chease_binary(config::EquilibriumConfig) # 1D array allocation zcpr, zcppr, zq, zdq, ztmf, ztp, zfb, zfbp, zpsi, zpsim = - [zeros(i == 1 || i == 10 ? npsi1-1 : npsi1) for i in 1:10] + [zeros(i == 1 || i == 10 ? npsi1 - 1 : npsi1) for i in 1:10] for arr in (zcpr, zcppr, zq, zdq, ztmf, ztp, zfb, zfbp, zpsi, zpsim) read(io, UInt32) @@ -275,7 +277,7 @@ function read_chease_ascii(config::EquilibriumConfig) data = Float64[] for line in lines[lines_range] for i in 0:4 - s = strip(line[(22*i+1):min(end, 22*(i+1))]) + s = strip(line[(22*i+1):min(end, 22 * (i + 1))]) if !isempty(s) push!(data, parse(Float64, s)) end @@ -405,7 +407,7 @@ function build_direct_from_ingest(config::EquilibriumConfig, ingest::DirectInges sq_in = cubic_interp(ingest.sq_xs, Series(ingest.sq_fs); extrap=ExtendExtrap()) psi_in = cubic_interp((ingest.psi_xs, ingest.psi_ys), ingest.psi_rz; extrap=ExtendExtrap()) return DirectRunInput(config, sq_in, psi_in, ingest.psi_xs, ingest.psi_ys, - ingest.rmin, ingest.rmax, ingest.zmin, ingest.zmax, ingest.psio, ingest.bt_sign, ingest) + ingest.rmin, ingest.rmax, ingest.zmin, ingest.zmax, ingest.psio, ingest.bt_sign, ingest.ip_sign, ingest) end """ @@ -432,15 +434,19 @@ Load an equilibrium from an IMAS data dictionary and return a `DirectRunInput`. The `dd.equilibrium.time_slice[]` is used (active time slice). Poloidal flux is converted from the IMAS COCOS convention (set by `config.imas_cocos`) to the internal COCOS 2 convention: + - `imas_cocos = 11` (default, IMAS standard): divide ψ by 2π - `imas_cocos = 2` (GPEC internal): no conversion ## Arguments -- `config`: `EquilibriumConfig` with `eq_type = "imas"` and `imas_cocos` set. -- `dd`: populated `IMASdd.dd` with `dd.equilibrium.time_slice[]` containing: - - `global_quantities.psi_axis`, `global_quantities.psi_boundary` - - `profiles_1d.psi`, `profiles_1d.f`, `profiles_1d.pressure`, `profiles_1d.q` - - `profiles_2d[1].grid.dim1` (R), `profiles_2d[1].grid.dim2` (Z), `profiles_2d[1].psi` + + - `config`: `EquilibriumConfig` with `eq_type = "imas"` and `imas_cocos` set. + + - `dd`: populated `IMASdd.dd` with `dd.equilibrium.time_slice[]` containing: + + + `global_quantities.psi_axis`, `global_quantities.psi_boundary` + + `profiles_1d.psi`, `profiles_1d.f`, `profiles_1d.pressure`, `profiles_1d.q` + + `profiles_2d[1].grid.dim1` (R), `profiles_2d[1].grid.dim2` (Z), `profiles_2d[1].psi` """ function read_imas(config::EquilibriumConfig, dd) @info "Processing IMAS equilibrium at global_time = $(dd.global_time) s" @@ -480,6 +486,10 @@ function read_imas(config::EquilibriumConfig, dd) # Capture toroidal-field sign from the boundary F value before abs() below. bt_sign = isempty(f_1d) ? 1 : Int(sign(f_1d[end])) bt_sign == 0 && (bt_sign = 1) + # Plasma-current sign from the IMAS global quantity (missing or zero → +1). + # The current is optional in IMAS, and reading an absent one throws rather than returning missing. + ip_imas = getproperty(eqt.global_quantities, :ip, 0.0) + ip_sign = ip_imas == 0 ? 1 : Int(sign(ip_imas)) nw = length(psi_1d) psi_norm_grid = range(0.0, 1.0; length=nw) @@ -528,7 +538,7 @@ function read_imas(config::EquilibriumConfig, dd) "\n Z ∈ [$(round(zmin; sigdigits=4)), $(round(zmax; sigdigits=4))] m" # Capture the raw arrays so an IMAS run can be replayed from gpec.h5 without the dd source. - ingest = DirectIngest(sq_xs, sq_fs_nodes, psi_in_xs, psi_in_ys, Matrix(psi_proc), rmin, rmax, zmin, zmax, psio, bt_sign) + ingest = DirectIngest(sq_xs, sq_fs_nodes, psi_in_xs, psi_in_ys, Matrix(psi_proc), rmin, rmax, zmin, zmax, psio, bt_sign, ip_sign) - return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, bt_sign, ingest) + return DirectRunInput(config, sq_in, psi_in, psi_in_xs, psi_in_ys, rmin, rmax, zmin, zmax, psio, bt_sign, ip_sign, ingest) end diff --git a/src/ForceFreeStates/Fourfit.jl b/src/ForceFreeStates/Fourfit.jl index 0f8095e3c..57cc28128 100644 --- a/src/ForceFreeStates/Fourfit.jl +++ b/src/ForceFreeStates/Fourfit.jl @@ -443,6 +443,7 @@ rational's resolution window, preserving the Δ′-stencil structure the equilib (`Equilibrium.RATIONAL_RES_RADIUS`). """ function core_capped_knots(xs::Vector{Float64}, rationals::Vector{Float64})::Vector{Int} + length(xs) < 3 && return collect(eachindex(xs)) cap_edge = 0.1 isempty(rationals) || (cap_edge = min(cap_edge, minimum(rationals) - Equilibrium.RATIONAL_RES_RADIUS)) in_rational_window(x) = any(abs(x - r) <= Equilibrium.RATIONAL_RES_RADIUS for r in rationals) diff --git a/src/ForceFreeStates/Kinetic.jl b/src/ForceFreeStates/Kinetic.jl index 0cd8e01dd..d148b12e9 100644 --- a/src/ForceFreeStates/Kinetic.jl +++ b/src/ForceFreeStates/Kinetic.jl @@ -1,3 +1,73 @@ + +# Knots across the envelope's [ψ_c, 2ψ_c] quintic-smoothstep transition, enough for a cubic spline +# to follow it without overshoot even on coarse decks. +const BAND_KNOTS = 9 + +""" + refine_grid_at_fbar_peaks(xs, kw, kt, evaluate, kin, equil, intr, psi_c; + ngrid=1000, relaxed_frac=0.01, target=3, max_add=24) → (xs, kw, kt) + +Insert kinetic evaluation knots across near-singular structure of F̄ that the grid does not +resolve. Scans cond(F̄) (the same operator `find_kinetic_singular_surfaces!` searches — Park & +Logan Eq. 70, so shifted and split resonances are included), takes peaks between +`relaxed_frac`·threshold and the singular threshold, measures each peak's FWHM, and adds knots +only where fewer than `target` knots lie inside it. New points respect `MIN_KNOT_SPACING`, stay +above the near-axis validity band, and are capped at `max_add`; each costs one kernel evaluation +and the existing values are reused. A well-resolved grid inserts nothing. +""" +function refine_grid_at_fbar_peaks(xs::Vector{Float64}, kw::Array{ComplexF64,3}, kt::Array{ComplexF64,3}, + evaluate::Function, kin::KineticMatrices, equil::Equilibrium.PlasmaEquilibrium, + intr::ForceFreeStatesInternal, psi_c::Float64; + ngrid::Int=1000, relaxed_frac::Float64=KINETIC_RELAXED_FRAC, target::Int=3, max_add::Int=24, + cond_threshold::Float64=KINETIC_SINGULAR_COND) + + lo, hi = xs[1], xs[end] + scan = collect(range(lo, hi; length=ngrid)) + hint = Ref(1) + cond_vals = [ + try + evaluate_fbar_condition(x, kin, equil, intr; hint=hint) + catch + Inf + end for x in scan + ] + + add = Float64[] + for i in 2:(ngrid-1) + c = cond_vals[i] + (c > cond_vals[i-1] && c > cond_vals[i+1] && relaxed_frac * cond_threshold < c <= cond_threshold) || continue + l = i + while l > 1 && cond_vals[l] > c / 2 + l -= 1 + end + r = i + while r < ngrid && cond_vals[r] > c / 2 + r += 1 + end + inside = count(x -> scan[l] <= x <= scan[r], xs) + inside >= target && continue + w = (scan[r] - scan[l]) / 3 + for x in (scan[i], scan[i] - w, scan[i] + w) + (lo < x < hi && x > 2 * psi_c) || continue + any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, xs) && continue + any(y -> abs(y - x) < Equilibrium.MIN_KNOT_SPACING, add) && continue + push!(add, x) + end + end + isempty(add) && return xs, kw, kt + length(add) > max_add && (add = sort(add)[1:max_add]) + + sort!(add) + @info "Kinetic grid: $(length(add)) knot(s) added across unresolved near-singular F̄ structure at " * + "ψ=$(round.(add; digits=4)) (cond peaks below the singular threshold)" + kw_new, kt_new = evaluate(add) + allxs = vcat(xs, add) + perm = sortperm(allxs) + kw_all = cat(kw, kw_new; dims=1)[perm, :, :] + kt_all = cat(kt, kt_new; dims=1)[perm, :, :] + return allxs[perm], kw_all, kt_all +end + """ build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; calculated_source=nothing) @@ -27,11 +97,17 @@ function build_kinetic_matrix_splines( mats::MatrixSplines, intr::ForceFreeStatesInternal, metric::MetricData; - calculated_source::Union{Nothing,Function}=nothing + calculated_source::Union{Nothing,Function}=nothing, + axis_validity_psi_c::Float64=0.0 ) xs = metric.xs mpsi = length(xs) + # Coarse equilibrium grids cannot represent env·(increment) across the validity band, and the + # spline overshoot can land on a rational; pin the band ends (C² there) and resolve the middle. + band_knots(lo, hi) = axis_validity_psi_c > 0 ? + [x for x in range(axis_validity_psi_c, 2 * axis_validity_psi_c; length=BAND_KNOTS) if lo < x < hi] : Float64[] + # Get raw kinetic matrices (scaling is baked into each source) if ctrl.kinetic_source == "fixed" kw_flat, kt_flat = fixed_kinetic_matrices(intr.mpert, intr.numpert_total, mpsi, ctrl.kinetic_factor, intr.mlow, mats, xs) @@ -42,7 +118,14 @@ function build_kinetic_matrix_splines( "calling build_kinetic_matrix_splines directly, or pass " * "`calculated_source=KineticForces.compute_calculated_kinetic_matrices` explicitly." ) - kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats) + band = band_knots(xs[1], xs[end]) + if isempty(band) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats) + else + xs = sort!(unique!(vcat(collect(xs), band))) + mpsi = length(xs) + kw_flat, kt_flat = calculated_source(ctrl, equil, intr, metric, mats; psis=xs) + end kw_flat .*= ctrl.kinetic_factor kt_flat .*= ctrl.kinetic_factor else @@ -50,7 +133,24 @@ function build_kinetic_matrix_splines( end # Pre-compute FKG derived matrices (corresponds to Fortran method=0) - return _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat) + mats = _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat; xs=xs) + + # The FKG splines now exist, so F̄ can be scanned: add knots only where near-singular structure + # (shifted/split kinetic resonances) falls in an interval that does not resolve it. + if ctrl.kinetic_source == "calculated" && calculated_source !== nothing && mats.kinetic !== nothing + xs2, kw_flat, kt_flat = refine_grid_at_fbar_peaks( + collect(xs), kw_flat, kt_flat, + psis -> begin + kwn, ktn = calculated_source(ctrl, equil, intr, metric, mats; psis=psis) + (kwn .* ctrl.kinetic_factor, ktn .* ctrl.kinetic_factor) + end, + mats.kinetic, equil, intr, axis_validity_psi_c) + if length(xs2) != length(xs) + xs = xs2 + mats = _compute_fkg_matrices(mats, equil, intr, metric, kw_flat, kt_flat; xs=xs) + end + end + return mats end """ @@ -73,9 +173,9 @@ function _compute_fkg_matrices( intr::ForceFreeStatesInternal, metric::MetricData, kw_flat::Array{ComplexF64,3}, - kt_flat::Array{ComplexF64,3} + kt_flat::Array{ComplexF64,3}; + xs::Vector{Float64}=metric.xs ) - xs = metric.xs mpsi = length(xs) np = intr.numpert_total mpert = intr.mpert diff --git a/src/ForceFreeStates/Surfaces/Asymptotics.jl b/src/ForceFreeStates/Surfaces/Asymptotics.jl index 26746c0ca..87cac444a 100644 --- a/src/ForceFreeStates/Surfaces/Asymptotics.jl +++ b/src/ForceFreeStates/Surfaces/Asymptotics.jl @@ -74,7 +74,7 @@ function compute_sing_asymptotics( # This is the parameter α but for all modes - α = 0 for non-resonant modes power[ipert_res] .= -alpha - power[ipert_res .+ intr.numpert_total] .= alpha + power[ipert_res.+intr.numpert_total] .= alpha # Zeroth-order non-resonant solutions for ipert in 1:intr.numpert_total @@ -111,7 +111,7 @@ function compute_sing_asymptotics( msg *= @sprintf(" m0mat(1,2)= %+.12e %+.12ei\n", real(m0mat[1, 2]), imag(m0mat[1, 2])) msg *= @sprintf(" m0mat(2,1)= %+.12e %+.12ei\n", real(m0mat[2, 1]), imag(m0mat[2, 1])) msg *= @sprintf(" m0mat(2,2)= %+.12e %+.12ei\n", real(m0mat[2, 2]), imag(m0mat[2, 2])) - di = m0mat[1, 1]*m0mat[2, 2] - m0mat[2, 1]*m0mat[1, 2] + di = m0mat[1, 1] * m0mat[2, 2] - m0mat[2, 1] * m0mat[1, 2] msg *= @sprintf(" di= %+.12e, alpha= %+.12e %+.12ei\n", real(di), real(alpha[1]), imag(alpha[1])) msg *= @sprintf(" psifac= %+.12e, r1=%d, ipert0=%d\n", singp.psifac, r1[1], ipert0) msg *= @sprintf(" vmat(ip,ip,2,0)= %+.8e %+.8ei\n", real(vmat[ipert0, ipert0, 2, 1]), imag(vmat[ipert0, ipert0, 2, 1])) @@ -579,7 +579,7 @@ function sing_get_ua(sing_asymp::SingAsymptotics, dpsi::Float64) # Restore powers (unshear v→u) — matches Fortran STRIDE sing_get_ua for i in eachindex(r1) - pfac = pfac_base ^ sing_asymp.alpha[i] # dpsi^α + pfac = pfac_base^sing_asymp.alpha[i] # dpsi^α ua[:, r2[2*i-1], :] ./= pfac # big solution column: /dpsi^α ua[:, r2[2*i], :] .*= pfac # small solution column: *dpsi^α ua[r1[i], :, 1] ./= sqrtfac # resonant row ξ: /√dpsi diff --git a/src/ForceFreeStates/Surfaces/Finding.jl b/src/ForceFreeStates/Surfaces/Finding.jl index 2c7bdac14..9edfe6e0d 100644 --- a/src/ForceFreeStates/Surfaces/Finding.jl +++ b/src/ForceFreeStates/Surfaces/Finding.jl @@ -125,9 +125,11 @@ function sing_lim!(intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, # strategy. Multi-n runs are not supported — the "outermost rational + dmlim/n" cutoff depends # on which n is used — and fall back to qhigh / psihigh truncation with a warning. if ctrl.set_psilim_via_dmlim && intr.nlow <= 0 - error("sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " * - "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " * - "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead.") + error( + "sing_lim!: set_psilim_via_dmlim = true requires a resolved toroidal range, but got intr.nlow=$(intr.nlow). " * + "Assign intr.nlow / intr.nhigh (from ctrl.nn_low / ctrl.nn_high) before calling sing_lim!, " * + "or set set_psilim_via_dmlim = false to truncate via qhigh / psihigh instead." + ) elseif ctrl.set_psilim_via_dmlim && intr.nlow != intr.nhigh @warn "set_psilim_via_dmlim = true is ignored for multi-n runs (nn_low=$(intr.nlow), nn_high=$(intr.nhigh)); falling back to qhigh / psihigh truncation." elseif ctrl.set_psilim_via_dmlim @@ -239,6 +241,11 @@ function evaluate_fbar_condition(psi::Float64, kin::KineticMatrices, equil::Equi return cond(fbar) end +# Kinetic F̄ counts as singular above this condition number; structure within +# KINETIC_RELAXED_FRAC of it is reported as near-singular (the shifted/split resonances). +const KINETIC_SINGULAR_COND = 1.0e8 +const KINETIC_RELAXED_FRAC = 0.01 + """ find_kinetic_singular_surfaces!(mats, equil, intr; ngrid=2000, cond_threshold=1e8) @@ -250,14 +257,16 @@ Mirrors the intent of Fortran `ksing_find` (`sing.f:1486-1616`) which finds zero det(F̄) via adaptive bisection. Here we use condition number peaks instead of determinant zeros for better numerical robustness and scale invariance. -Algorithm: - - 1. Evaluate cond(F̄) on a dense ψ grid - 2. Find local maxima (peaks where gradient changes from + to -) - 3. Refine each peak with golden-section minimization of -cond - 4. Filter by threshold and resonance condition +Evaluates cond(F̄) on a dense ψ grid, finds its local maxima, refines each with golden-section +minimization of -cond, then filters by threshold and resonance condition. """ -function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium.PlasmaEquilibrium, intr::ForceFreeStatesInternal; ngrid::Int=2000, cond_threshold::Float64=1e8) +function find_kinetic_singular_surfaces!( + mats::MatrixSplines, + equil::Equilibrium.PlasmaEquilibrium, + intr::ForceFreeStatesInternal; + ngrid::Int=2000, + cond_threshold::Float64=KINETIC_SINGULAR_COND +) kin = mats.kinetic kin === nothing && error("find_kinetic_singular_surfaces! requires a kinetic fit; call build_kinetic_matrix_splines first") psilow = equil.profiles.xs[1] @@ -282,11 +291,17 @@ function find_kinetic_singular_surfaces!(mats::MatrixSplines, equil::Equilibrium intr.kinsing_scan_threshold = cond_threshold # Find local maxima of cond(F̄): points where cond increases then decreases - peak_indices = Int[] - for i in 2:(ngrid-1) - if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1] && cond_vals[i] > cond_threshold - push!(peak_indices, i) - end + local_maxima = [i for i in 2:(ngrid-1) if cond_vals[i] > cond_vals[i-1] && cond_vals[i] > cond_vals[i+1]] + peak_indices = filter(i -> cond_vals[i] > cond_threshold, local_maxima) + + # Peaks below threshold mark near-singular structure (the shifted/split resonances of Park & + # Logan Eq. 70); report the strongest few so sharp kinetic structure is visible rather than silent. + subthreshold = filter(i -> KINETIC_RELAXED_FRAC * cond_threshold < cond_vals[i] <= cond_threshold, local_maxima) + if !isempty(subthreshold) + top = sort(subthreshold; by=i -> cond_vals[i], rev=true)[1:min(3, length(subthreshold))] + @info "Kinetic F̄ near-singular structure below the singular threshold at " * + join(["ψ=$(round(psi_grid[i]; digits=4)) (cond=$(round(cond_vals[i]; sigdigits=3)))" for i in top], ", ") * + " — full scan in SingularSurfaces/Kinetic/scan_cond; check the ψ grid resolves these if results look grid-sensitive" end # Refine each peak to find the precise ψ location diff --git a/src/ForcingTerms/CoilFourier.jl b/src/ForcingTerms/CoilFourier.jl index 62996bf59..9b2c5b22c 100644 --- a/src/ForcingTerms/CoilFourier.jl +++ b/src/ForcingTerms/CoilFourier.jl @@ -100,14 +100,15 @@ via periodic cubic splines on the resulting R(θ_norm), Z(θ_norm) data. The toroidal grid direction follows the Fortran GPEC convention: `phi_j = -helicity × 2π × j/nzeta` where `helicity = sign(Bt) × sign(Ip)`. -This is derived from `equil.params.bt_sign` and `equil.params.crnt`. +This is derived from `equil.params.bt_sign` and `equil.params.ip_sign` (the signs the equilibrium +file states; the computed current `crnt` is always positive and carries no direction). For DIII-D (Bt < 0, Ip > 0 → helicity = -1): phi increases with j (standard direction). For positive-helicity machines (Bt > 0, Ip > 0 → helicity = +1): phi decreases with j. """ function sample_boundary_grid(equil::Equilibrium.PlasmaEquilibrium, mtheta::Int, nzeta::Int; psi::Float64=equil.rzphi_xs[end]) # Build uniform theta grid (same convention as equil.rzphi_ys, but potentially finer) - theta_grid = range(0; length=mtheta, step=1.0/mtheta) + theta_grid = range(0; length=mtheta, step=1.0 / mtheta) R_arr = zeros(mtheta) Z_arr = zeros(mtheta) @@ -136,10 +137,8 @@ function sample_boundary_grid(equil::Equilibrium.PlasmaEquilibrium, mtheta::Int, # Helicity sets the direction of the toroidal angle grid to match Fortran convention: # phi_j = -helicity × 2π × j/nzeta, helicity = sign(Bt) × sign(Ip) - bt_sign = !isnothing(equil.params.bt_sign) ? equil.params.bt_sign : 1 - ip_sign = !isnothing(equil.params.crnt) ? Int(sign(equil.params.crnt)) : 1 - helicity = bt_sign * ip_sign - phi_grid = collect(range(0; length=nzeta, step=(-helicity * 2π/nzeta))) + helicity = equil.params.bt_sign * equil.params.ip_sign + phi_grid = collect(range(0; length=nzeta, step=(-helicity * 2π / nzeta))) # Toroidal angle offset ν(ψ, θ_SFL): in SFL coordinates the physical toroidal angle at # grid point (θ_SFL, ζ_SFL) is φ_phys = -helicity*(2π*ζ_SFL + ν(ψ,θ_SFL)). diff --git a/src/GeneralizedPerturbedEquilibrium.jl b/src/GeneralizedPerturbedEquilibrium.jl index f2ed4ef65..df8a0678e 100755 --- a/src/GeneralizedPerturbedEquilibrium.jl +++ b/src/GeneralizedPerturbedEquilibrium.jl @@ -85,8 +85,8 @@ using .Equilibrium: PlasmaEquilibrium using .ForcingTerms: RMPField const _DEPRECATED_FFS_KEYS = ("mer_flag", "force_wv_symmetry", "ode_flag", "cyl_flag", "mat_flag", "reform_eq_with_psilim", - "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi", - "gal_flag") + "use_riccati", "use_parallel", "parallel_threads", "populate_dense_xi", + "gal_flag") const _DEPRECATED_EQUIL_KEYS = ("power_bp", "power_b", "power_r", "power_rc") # Drop deprecated keys from a parsed gpec.toml section so legacy files keep parsing @@ -210,7 +210,8 @@ function main_from_inputs( equil = Equilibrium.setup_equilibrium(eq_config, additional_input) kf_ctrl, kinetic_profiles, kf_species = load_kinetic_context(inputs, intr, ctrl, equil) - equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) + equil = maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; + kf_ctrl=kf_ctrl, kf_species=kf_species) @info "Equilibrium construction completed in $(@sprintf("%.3f", time() - equil_start)) s" @@ -256,6 +257,15 @@ function main_from_inputs( @info "Results written to $(ctrl.HDF5_filename)" end + # Validity describes the kinetic model, not the torque contraction, so it is written wherever + # kinetic profiles were loaded — including FFS-only runs that never reach run_kinetic_forces. + if ctrl.write_outputs_to_HDF5 && kinetic_profiles !== nothing + h5open(joinpath(ffs_result.dir_path, ctrl.HDF5_filename), "cw") do h5file + KineticForces.write_validity!(h5file, kf_ctrl, kf_species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) + end + end + @info "Force-Free States completed in $(@sprintf("%.3f", time() - ffs_start)) s" # Early exit if user only requested force-free states (SLAYER still runs). @@ -379,7 +389,8 @@ function load_kinetic_context( end """ - maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles) -> equil + maybe_reform_equilibrium(equil, eq_config, additional_input, intr, ctrl, kinetic_profiles; + kf_ctrl=nothing, kf_species=nothing) -> equil Two-pass auto grid: measure the pass-1 equilibrium's curvature (profiles, geometry, kinetic profiles), pin knots on rational surfaces, and re-form on the refined grid from the in-memory @@ -391,16 +402,24 @@ function maybe_reform_equilibrium( additional_input, intr::ForceFreeStatesInternal, ctrl::ForceFreeStatesControl, - kinetic_profiles + kinetic_profiles; + kf_ctrl::Union{Nothing,KineticForces.KineticForcesControl}=nothing, + kf_species=nothing ) Equilibrium.wants_two_pass(eq_config) || return equil mandatory = ForceFreeStates.rational_psi_nodes(equil; nlow=intr.nlow, nhigh=intr.nhigh) + # Kinetic-resonance nodes are located here rather than at the call site so the validity + # boundary that filters them is cleared against this pass's rationals (intr.sing is empty yet). + pinned = kf_ctrl === nothing ? Float64[] : + KineticForces.resonance_grid_nodes(ctrl, kf_ctrl, kinetic_profiles, kf_species, equil, intr; rationals=mandatory) # Smallest |n| in the run sets the widest matching half-stencil dpsi = singfac_min/(n_min·|q′|), # so the rational-surface brackets clear a zone large enough for every mode. n_min = minimum(abs(n) for n in intr.nlow:intr.nhigh if n != 0) + isempty(pinned) || + @info "Pinning $(length(pinned)) kinetic-resonance surfaces into the ψ grid: $(round.(sort(pinned); digits=3))" psi_nodes = Equilibrium.refined_psi_grid(equil; - tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, + tau=eq_config.psi_accuracy, kin=kinetic_profiles, mandatory=mandatory, pinned=pinned, singfac_min=ctrl.singfac_min, n_min=n_min) rerun_input = if additional_input !== nothing # Analytic *Config, IMAS dd, or prebuilt RunInput — all re-formable. The IMAS @@ -579,12 +598,18 @@ function prepare_force_free_states!( # Inject the KineticForces callback so the "calculated" source can # invoke compute_calculated_kinetic_matrices without ForceFreeStates # importing KineticForces (which would invert the load order). - 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, species=species) + # One boundary for the whole run: the widest-orbit species, moved clear of any rational + # surface whose window the transition band would otherwise cut through. + axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, equil, + Float64[sng.psifac for sng in intr.sing]) + calculated_cb = + (c, e, i, m, f; psis::Vector{Float64}=Float64[]) -> + KineticForces.compute_calculated_kinetic_matrices( + c, e, i, m, f; + kf_ctrl=kf_ctrl, kinetic_profiles=kinetic_profiles, species=species, psis=psis, + axis_psi_c=axis_psi_c) mats = build_kinetic_matrix_splines(ctrl, equil, mats, intr, metric; - calculated_source=calculated_cb) + calculated_source=calculated_cb, axis_validity_psi_c=axis_psi_c) # Find kinetically-displaced singular surfaces (zeros of det(F̄)) for ODE crossings. # Matches Fortran ksing_find (sing.f:1486-1616). singfac_min > 0 gates crossings; @@ -740,10 +765,10 @@ runs are TOML-driven this cycle: `kinetic_factor > 0` needs the `[KineticForces] and errors here. ```julia -eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") +eq = PlasmaEquilibrium("input.geqdsk"; jac_type="hamada") prob = EulerLagrangeProblem(eq; nn=1, delta_mlow=8, delta_mhigh=8, vac_flag=true) -ffs = solve(prob, Riccati()) -ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # equivalent one-line form +ffs = solve(prob, Riccati()) +ffs = solve(eq, Riccati(); nn=1, vac_flag=true) # equivalent one-line form ``` """ function solve(prob::EulerLagrangeProblem, alg::ForceFreeStates.AbstractIntegrator) @@ -871,7 +896,7 @@ function perturbed_equilibrium( kwargs... ) ctrl = ffs.control - pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; kwargs...) + pe_ctrl = PerturbedEquilibrium.PerturbedEquilibriumControl(; PerturbedEquilibrium.kinetic_regularization_kwargs(ffs, kwargs)...) pe_intr = PerturbedEquilibrium.PerturbedEquilibriumInternal(; dir_path=ffs.dir_path) # Inner-layer penetrated resonant field; zeros under ideal closure. @@ -932,6 +957,10 @@ function run_kinetic_forces( @info "\n KineticForces\n$_SECTION" kf_start = time() + # One boundary for the whole run, shared with the EL kinetic matrices and the Validity output. + axis_psi_c = KineticForces.axis_validity_boundary(kf_ctrl, species, kinetic_profiles, result.equil, + Float64[sng.psifac for sng in result.surfaces]) + # Standalone NTV torque diagnostics need a PE state (they contract kinetic operators # against ξ). The self-consistent kinetic_source="calculated" path produces none — skip. if pe_state === nothing @@ -943,7 +972,8 @@ function run_kinetic_forces( if species === nothing kf_state = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles) + KineticForces.compute_torque_all_methods!(kf_state, kf_intr, kf_ctrl, result.equil, kinetic_profiles; + axis_psi_c=axis_psi_c) if kf_ctrl.write_outputs_to_HDF5 h5open(joinpath(result.dir_path, kf_ctrl.HDF5_filename), "cw") do h5file KineticForces.write_to_hdf5!(h5file, kf_state; @@ -960,7 +990,8 @@ function run_kinetic_forces( (f => getfield(kf_ctrl, f) for f in fieldnames(KineticForces.KineticForcesControl))..., zi=sp.z, mi=sp.m, electron=sp.electron, ion_species=KineticForces.IonSpecies[]) st = KineticForces.KineticForcesState() - KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles) + KineticForces.compute_torque_all_methods!(st, kf_intr, sctrl, result.equil, sp.profiles; + axis_psi_c=axis_psi_c) push!(states, st) push!(labels, sp.label) end diff --git a/src/HDF5Schema.jl b/src/HDF5Schema.jl index 274f54f14..2f2f6f175 100644 --- a/src/HDF5Schema.jl +++ b/src/HDF5Schema.jl @@ -44,6 +44,7 @@ const EQUIL_H5_NAMES = Dict( :li2 => "l_i_2", :li3 => "l_i_3", :bt_sign => "B_T_sign", + :ip_sign => "I_p_sign", :psio => "psi_total" ) const EQUIL_H5_SKIP = Set([:psi0, :psi_axis, :psi_axis_norm, :zsep, :verbose, :diagnose_src, :diagnose_maxima]) @@ -98,6 +99,7 @@ const MAIN_H5_ANNOTATIONS = [ "Equilibrium/l_i_3" => (; long_name="internal inductance (definition 3)"), "Equilibrium/volume" => (; long_name="plasma volume", units="m^3"), "Equilibrium/B_T_sign" => (; long_name="sign of the toroidal field"), + "Equilibrium/I_p_sign" => (; long_name="sign of the plasma current as stated by the equilibrium file"), "Equilibrium/psi_norm" => (; long_name="normalized poloidal flux at the magnetic axis (0 by definition of ψ_N)"), "Equilibrium/psi_boundary" => (; long_name="poloidal flux at the plasma boundary in the internal normalized convention (1 by construction, not a Wb/rad datum)"), "Equilibrium/psi_boundary_norm" => (; long_name="normalized poloidal flux at the plasma boundary (1 by definition of ψ_N)"), @@ -245,7 +247,7 @@ const MAIN_H5_ANNOTATIONS = [ "SurfaceGeometries/Plasma/z" => (; long_name="Cartesian z of plasma-surface point cloud", units="m"), "SurfaceGeometries/Wall/x" => (; long_name="Cartesian x of wall point cloud", units="m"), "SurfaceGeometries/Wall/y" => (; long_name="Cartesian y of wall point cloud", units="m"), - "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m"), + "SurfaceGeometries/Wall/z" => (; long_name="Cartesian z of wall point cloud", units="m") ] # Euler-Lagrange operator matrices: same wording per letter, Ideal/ and Kinetic/ variants. @@ -258,7 +260,7 @@ const _ELM_IDEAL_LETTERS = [ ("H", "Euler-Lagrange primitive coefficient matrix H"), ("F", "Euler-Lagrange derived coefficient matrix F"), ("K", "Euler-Lagrange derived coefficient matrix K"), - ("G", "Euler-Lagrange derived coefficient matrix G"), + ("G", "Euler-Lagrange derived coefficient matrix G") ] # The kinetic branch overwrites only A, B, C, K, G and adds f0; D, E, H, F are shared # unchanged from the ideal set and are not re-emitted. @@ -268,14 +270,19 @@ const _ELM_KINETIC_LETTERS = [ ("C", "Euler-Lagrange primitive coefficient matrix C"), ("K", "Euler-Lagrange derived coefficient matrix K"), ("G", "Euler-Lagrange derived coefficient matrix G"), - ("f0", "raw kinetic component matrix f0"), + ("f0", "raw kinetic component matrix f0") ] const ELM_H5_ANNOTATIONS = vcat( ["ForceFreeStates/EulerLagrangeMatrices/psi" => (; long_name="normalized poloidal flux ψ_N grid of the operator matrices", scale="psi")], - ["ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => - (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS], - ["ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => - (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_KINETIC_LETTERS] + [ + "ForceFreeStates/EulerLagrangeMatrices/Ideal/$l" => + (; long_name="ideal " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for (l, d) in _ELM_IDEAL_LETTERS + ], + [ + "ForceFreeStates/EulerLagrangeMatrices/Kinetic/$l" => + (; long_name="kinetic-modified " * d, dims=("psi", "mode_row", "mode_col"), attach=(1 => "ForceFreeStates/EulerLagrangeMatrices/psi",)) for + (l, d) in _ELM_KINETIC_LETTERS + ] ) """ diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index 43d1fc214..ddbde66c1 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -25,7 +25,7 @@ import .GGJ: delta_convergence, solution_profile, asymptotic_profile, q4_surface import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths -import .SLAYER: surface_minor_radius, surface_da_dpsi, build_slayer_inputs +import .SLAYER: surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner, solve_inner_profile export GGJ, GGJModel, GGJParameters @@ -37,6 +37,6 @@ export delta_convergence, solution_profile, asymptotic_profile, q4_surface_bench export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths -export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs end # module InnerLayer diff --git a/src/InnerLayer/SLAYER/LayerInputs.jl b/src/InnerLayer/SLAYER/LayerInputs.jl index 96177903d..252e81ae9 100644 --- a/src/InnerLayer/SLAYER/LayerInputs.jl +++ b/src/InnerLayer/SLAYER/LayerInputs.jl @@ -9,14 +9,14 @@ # Geometry extraction: # - Minor radius at the outboard midplane (θ = 0) via # `equil.rzphi_rsquared((ψ, 0.0))`. -# - `da/dψ` via central finite difference on the same bicubic. +# - `da/dψ` from the interpolant's own analytic ψ-derivative. # - r-based magnetic shear via `r_based_shear(rs, q, q1, da/dψ)` (defined # in LayerParameters.jl). using ..Utilities: KineticProfiles using ...Utilities.NeoclassicalResistivity: NeoResistivityModel, SpitzerModel, coulomb_log_e, nu_star_e -using FastInterpolations: DerivOp +using FastInterpolations: DerivOp, integrate, cubic_interp, cumulative_integrate, ExtendExtrap """ surface_minor_radius(equil, psi; theta=0.0) -> Float64 @@ -32,33 +32,126 @@ function surface_minor_radius(equil, psi::Real; theta::Real=0.0) end """ - surface_da_dpsi(equil, psi; theta=0.0, h=1e-5) -> Float64 + surface_da_dpsi(equil, psi; theta=0.0) -> Float64 -Central finite-difference approximation of `d(minor radius)/dψ` at `psi`. -Falls back to one-sided differences near the flux-coordinate boundaries -(0 or 1). +Analytic ψ-derivative of the minor radius at `psi` and poloidal angle +`theta`, taken from the `rzphi_rsquared` interpolant's own ψ-derivative as +`da/dψ = (∂r²/∂ψ) / (2a)`. Valid wherever the interpolant is, including +under extrapolation past the ψ grid. Diverges at the magnetic axis, where +`a ~ √ψ`; callers evaluating near `ψ = 0` must check `isfinite`. """ -function surface_da_dpsi(equil, psi::Real; theta::Real=0.0, h::Real=1e-5) - psi_f = Float64(psi) - # Clamp to safe sampling range within (0, 1) - eps_edge = 10 * h - lo = psi_f - h - hi = psi_f + h - if lo < eps_edge - # one-sided forward - a0 = surface_minor_radius(equil, max(psi_f, eps_edge); theta=theta) - a1 = surface_minor_radius(equil, max(psi_f, eps_edge) + h; theta=theta) - return (a1 - a0) / h - elseif hi > 1.0 - eps_edge - # one-sided backward - a0 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge) - h; theta=theta) - a1 = surface_minor_radius(equil, min(psi_f, 1.0 - eps_edge); theta=theta) - return (a1 - a0) / h +function surface_da_dpsi(equil, psi::Real; theta::Real=0.0) + return _da_dpsi_at_theta(equil, Float64(psi), Float64(theta)) +end + +# d(√r²)/dψ at one (ψ, θ) from the interpolant's own ψ-derivative. Shared by all radial-label +# conventions so none carries its own stencil. +@inline function _da_dpsi_at_theta(equil, psi::Float64, theta::Float64) + r_sq = equil.rzphi_rsquared((psi, theta)) + a = sqrt(max(r_sq, 0.0)) + a > 0 || return Inf # magnetic axis: a ~ √ψ, so da/dψ genuinely diverges + return equil.rzphi_rsquared((psi, theta); deriv=DerivOp(1, 0)) / (2a) +end + +""" + radial_label(equil; rs_method=:midplane, theta=0.0) -> (r_at, dr_dpsi_at) + +Build the pair of closures `r_at(ψ)` and `dr_dpsi_at(ψ)` defining one radial +label for the layer stack. Both closures must come from the same label, +because the r-based shear `(r/q)dq/dr`, `τ_R = μ₀r²/η`, `τ_E = r²/χ`, the +`d_β/r` normalization, the Δ' reference-length factor `k_ref = r_s/(da/dψ)`, +and any metre-to-ψ width conversion all have to live in one coordinate; a +mismatched `r` and `dr/dψ` silently corrupts every one of them. All +derivatives are analytic — no label carries a finite-difference stencil. + +# Labels + + - `:midplane` -- outboard-midplane chord from the magnetic axis at `theta` + (historical default), natural for comparison with midplane diagnostics. + - `:halfwidth` -- midplane half-chord, the mean of the outboard and inboard + chords at `θ = 0` and `θ = 0.5`; shift-free. Coincides with the flux label + on circular equilibria but is its own convention on shaped ones. + - `:fsa` -- θ-mean surface radius, a 128-point midpoint mean of the local + minor radius; the closest geometric approximation to `:flux` at interior + surfaces of shaped equilibria. + - `:volume` -- cylinder-equivalent label `√(V(ψ)/(2π²R₀))`, the + Rutherford-literature convention. + - `:flux` -- toroidal-flux label. Fitzpatrick, Nucl. Fusion (2025), + Eq. 30: `dψ_p/dr = B₀ r g/q` integrates to `ψ_t = B₀r²/2`, so + `r = √(2ψ_t/B₀)` with `ψ_t = psio·∫₀^ψ (q/g) dψ′` and `g = F/(B₀R₀)`, + `F = R·B_tor`. Note `profiles.F_spline` stores `2πF`, so the + implementation divides it by 2π to form `g`. + Defined from flux alone, it carries no circular-cross-section assumption, + and its derivative `dr/dψ ∝ q` grows toward a separatrix where the + geometric labels' `da/dψ` collapses. + +On shaped equilibria the labels agree at low-q surfaces and diverge strongly +near the edge, where the slab-layer matching is label-ambiguous regardless of +choice. The label is selected programmatically; it is not exposed via TOML. +""" +function radial_label(equil; rs_method::Symbol=:midplane, theta::Real=0.0) + theta_f = Float64(theta) + + _flux_r, _flux_dr = if rs_method === :flux + b0f = Float64(equil.params.b0) + R0f = Float64(equil.ro) + psiof = Float64(equil.psio) + xs_f = collect(Float64, equil.profiles.xs) + # Carry g = F/(B0 R0) rather than assuming g = 1: r² = 2∫(q/g)dψ_p/B0. + _g_at(x) = Float64(equil.profiles.F_spline(x)) / (2π * b0f * R0f) + qg = [Float64(equil.profiles.q_spline(x)) / _g_at(x) for x in xs_f] + Phi = collect(Float64, cumulative_integrate(cubic_interp(xs_f, qg))) + r_knots = sqrt.(max.(2 .* psiof .* Phi ./ b0f, 0.0)) + rspl = cubic_interp(xs_f, r_knots; extrap=ExtendExtrap()) + (ψ -> Float64(rspl(Float64(ψ))), + ψ -> psiof * Float64(equil.profiles.q_spline(Float64(ψ))) / + (b0f * _g_at(Float64(ψ)) * max(Float64(rspl(Float64(ψ))), eps()))) else - a_plus = surface_minor_radius(equil, psi_f + h; theta=theta) - a_minus = surface_minor_radius(equil, psi_f - h; theta=theta) - return (a_plus - a_minus) / (2h) + (nothing, nothing) end + + _a_at(ψ, θ) = sqrt(max(equil.rzphi_rsquared((Float64(ψ), Float64(θ))), 0.0)) + + _rs_at(ψ) = + if rs_method === :fsa + N = 128 + s = 0.0 + @inbounds for k in 1:N + s += _a_at(ψ, (k - 0.5) / N) + end + s / N + elseif rs_method === :halfwidth + 0.5 * (_a_at(ψ, 0.0) + _a_at(ψ, 0.5)) + elseif rs_method === :volume + # Lower bound just off the axis, where dV/dψ is extrapolated; the omitted sliver + # of core volume is negligible for an edge-surface label. + V = integrate(equil.profiles.dVdpsi_spline, 1e-4, Float64(ψ)) + sqrt(max(V, 0.0) / (2π^2 * equil.ro)) + elseif rs_method === :flux + _flux_r(ψ) + else + surface_minor_radius(equil, ψ; theta=theta_f) + end + + _da_dpsi_at(ψ) = + if rs_method === :fsa + N = 128 + s = 0.0 + @inbounds for k in 1:N + s += _da_dpsi_at_theta(equil, Float64(ψ), (k - 0.5) / N) + end + s / N + elseif rs_method === :halfwidth + 0.5 * (_da_dpsi_at_theta(equil, Float64(ψ), 0.0) + _da_dpsi_at_theta(equil, Float64(ψ), 0.5)) + elseif rs_method === :volume + Float64(equil.profiles.dVdpsi_spline(ψ)) / (4π^2 * equil.ro * max(_rs_at(ψ), eps())) + elseif rs_method === :flux + _flux_dr(ψ) + else + _da_dpsi_at_theta(equil, Float64(ψ), theta_f) + end + + return (_rs_at, _da_dpsi_at) end """ @@ -114,6 +207,12 @@ profiles, without an intermediate file round-trip. a prescribed value. (For `dc_type=:rfitzp` and `:lar`, dgeo_val is not consulted.) - `dc_type` -- `:none` (default), `:lar`, `:rfitzp`, or `:toroidal`. + - `rs_method` -- radial label defining `r_s` for the whole layer stack: + `:midplane` (default), `:halfwidth`, `:fsa`, `:volume`, or `:flux`. See + [`radial_label`](@ref) for the definitions. S, the r-based shear, W_d, + and the Δ' reference-length factor `k_ref` all follow the choice + together, so every option is self-consistent. Not exposed via TOML — + programmatic use only. - `theta` -- poloidal angle at which to measure minor radius (default `0.0`, outboard midplane). - `resistivity_model` -- `SauterNeoModel()` (default), `RedlNeoModel()`, @@ -156,54 +255,22 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; Float64(bt(ψ)) end - # Minor-radius extractor: `:midplane` = outboard-midplane chord - # (original behavior); `:fsa` = θ-mean of √rzphi_rsquared, the - # flux-surface-averaged minor radius. - _rs_at(ψ) = - if rs_method === :fsa - integrand(θ) = sqrt(equil.rzphi_rsquared((Float64(ψ), Float64(θ)))) - N = 128 - s = 0.0 - @inbounds for k in 1:N - s += integrand((k - 0.5) / N) - end - s / N - else - surface_minor_radius(equil, ψ; theta=theta) - end - _da_dpsi_at(ψ) = - if rs_method === :fsa - # central finite difference on _rs_at - h = 1e-5 - lo = ψ - h - hi = ψ + h - eps_edge = 10h - if lo < eps_edge - (_rs_at(max(ψ, eps_edge) + h) - _rs_at(max(ψ, eps_edge))) / h - elseif hi > 1.0 - eps_edge - (_rs_at(min(ψ, 1.0 - eps_edge)) - _rs_at(min(ψ, 1.0 - eps_edge) - h)) / h - else - (_rs_at(ψ + h) - _rs_at(ψ - h)) / (2h) - end - else - surface_da_dpsi(equil, ψ; theta=theta) - end + _rs_at, _da_dpsi_at = radial_label(equil; rs_method=rs_method, theta=theta) - # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline - # derivatives. When `compute_omega_star=true` we override any ω_*e/ω_*i - # carried in `profiles`. Main-ion density is - # taken equal to the electron density (quasi-neutrality, matching the - # staging step). + # Per-surface ω_*e, ω_*i (diamagnetic frequencies) from spline derivatives. + # In flux coordinates ω_* = n·(dp/dψ)/(e·n_e), so the surface's toroidal mode + # number enters here. Main-ion density is taken equal to the electron density + # (quasi-neutrality, matching the staging step). chi1 = 2π * equil.psio - _omega_star_at(ψ) = begin + _omega_star_at(ψ, n_tor) = begin n_e = Float64(profiles.n_e(ψ)) dn_e = Float64(profiles.n_e(ψ; deriv=DerivOp(1))) T_e = Float64(profiles.T_e(ψ)) dT_e = Float64(profiles.T_e(ψ; deriv=DerivOp(1))) T_i = Float64(profiles.T_i(ψ)) dT_i = Float64(profiles.T_i(ψ; deriv=DerivOp(1))) - ω_star_e = (2π / chi1) * (T_e * dn_e / n_e + dT_e) - ω_star_i = -(2π / (Float64(z_i) * chi1)) * (T_i * dn_e / n_e + dT_i) + ω_star_e = n_tor * (2π / chi1) * (T_e * dn_e / n_e + dT_e) + ω_star_i = -n_tor * (2π / (Float64(z_i) * chi1)) * (T_i * dn_e / n_e + dT_i) return (ω_star_e, ω_star_i) end @@ -217,20 +284,17 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; da_dpsi = _da_dpsi_at(psi) sval_r = r_based_shear(rs, q, q1, da_dpsi) - prof = profiles(psi) - # Override ω_*e, ω_*i with spline-derivative values when requested. - ω_e_use, ω_i_use = if compute_omega_star - _omega_star_at(psi) - else - (prof.omega_e, prof.omega_i) - end - # Resonant (m, n): take the first element of the mode-number vectors. # Parallel-FM `sing.m`/`sing.n` hold exactly one entry each; ideal # DCON may hold multiple — we pick the first and document the choice. m_res = sing.m[1] n_res = sing.n[1] + prof = profiles(psi) + # Take ω_*e, ω_*i from the spline derivatives, or from `profiles` when the caller + # supplies them directly. `run_slayer` supplies zeros, so the latter is a library path. + ω_e_use, ω_i_use = compute_omega_star ? _omega_star_at(psi, n_res) : (prof.omega_e, prof.omega_i) + # Pull geometric trapped-fraction inputs from ResistGeometry when # available (populated by ForceFreeStates.resist_eval_all!); else # fall back to nothing and let slayer_parameters compute them from @@ -289,6 +353,26 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; _eval(dgeo_val, psi) end + # Reference-length conversion inputs for the outer Δ': K = r_s·(dψ_N/dr)|_s and + # α = √(−D_I) (Glasser-Greene-Johnson 1975 Eq. 48), with α clamped to 0 on Mercier-unstable + # surfaces (the factor turns complex there) and K = 1 whenever da/dψ is not a usable + # positive number (zero/non-finite/negative would corrupt the Δ' diagonal). + k_ref_k = if isfinite(da_dpsi) && da_dpsi > 0.0 + rs / da_dpsi + else + @warn("build_slayer_inputs: da/dψ = $da_dpsi at ψ = $psi is not usable; leaving " * + "Δ' unconverted (k_ref = 1) at this surface.", maxlog=3) + 1.0 + end + alpha_k = if rg === nothing + @warn("build_slayer_inputs: sing.restype not populated; using the " * + "slab Mercier exponent α = 1/2 for the Δ' reference-length " * + "conversion at all such surfaces.", maxlog=1) + 0.5 + else + sqrt(max(-(rg.E + rg.F + rg.H - 0.25), 0.0)) + end + out[k] = slayer_parameters(; n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i, omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use, @@ -304,7 +388,9 @@ function build_slayer_inputs(equil, sings, profiles::KineticProfiles; f_trap=f_trap_kw, nu_e_star=nu_e_star_kw, R_major_eff=R_major_eff, - lnLambda_form=lnLambda_form + lnLambda_form=lnLambda_form, + k_ref=k_ref_k, + alpha_mercier=alpha_k ) end return out diff --git a/src/InnerLayer/SLAYER/LayerParameters.jl b/src/InnerLayer/SLAYER/LayerParameters.jl index 9718ad666..ae875dea4 100644 --- a/src/InnerLayer/SLAYER/LayerParameters.jl +++ b/src/InnerLayer/SLAYER/LayerParameters.jl @@ -27,7 +27,7 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and | `tau` | T_i / T_e | | `lu` | Lundquist number S = τ_R / τ_H | | `c_beta` | Compressibility √(β_local / (1 + β_local)) | -| `D_norm` | (d_β/r_s) · S^(1/3) · √(τ/(1+τ)) (Fitzpatrick normalized scale) | +| `D_norm` | (d_β/r_s) · S^(1/3) · √ι_e (Fitzpatrick normalized scale) | | `P_perp` | Perpendicular Prandtl number τ_R / τ_⊥ | | `P_tor` | Toroidal-direction Prandtl number τ_R / τ_‖tor | | `Q_e` | Normalized electron diamagnetic: −tauk · ω_*e | @@ -46,6 +46,15 @@ de-normalization. The parametrization uses `P_perp`, `P_tor`, and | `d_beta` | Beta-weighted ion length scale c_β · d_i [m] | | `dc_tmp` | Critical-Δ offset from chi_parallel matching | | `dc_type` | Selector for `dc_tmp` formula | +| `k_ref` | Reference-length ratio K = r_s · (dψ_N/dr) at this surface (1 = no Δ' conversion) | +| `alpha_mercier` | Mercier Frobenius exponent α = √(−D_I) (Glasser-Greene-Johnson 1975 Eq. 48) governing the Δ' reference-length conversion (1/2 = slab/cylindrical value) | + +`k_ref` and `alpha_mercier` feed the ψ_N → r_s reference-length conversion of +the outer Δ' matrix (see `delta_prime_to_rs_reference` in the Tearing +runner): the slab layer, its `dc_tmp` critical-Δ, and the `S^(1/3)` Δ(Q) +scale are all referenced to unit `x̂ = (r−r_s)/r_s`, while the outer BVP Δ' +is referenced to unit Δψ_N. Hand-built parameters default to `k_ref = 1`, +which makes the conversion the identity. The complex normalized growth rate `Q = ω + iγ` is **not** stored here; it is passed as a separate argument to `solve_inner`. @@ -85,6 +94,10 @@ Base.@kwdef struct SLAYERParameters <: InnerLayerParameters # Critical-Δ offset dc_tmp::Float64 = 0.0 dc_type::Symbol = :none + + # Reference-length conversion inputs for the outer Δ' (ψ_N → r_s-based x̂) + k_ref::Float64 = 1.0 + alpha_mercier::Float64 = 0.5 end # Allowed dc_type values for the critical-Δ offset. `:none` is the default @@ -202,6 +215,12 @@ parametrization (P_perp/P_tor/D_norm; the older magnetic/electron Prandtl - `dr_val`, `dgeo_val` -- inputs for the critical-Δ formula - `dc_type` -- one of `:none`, `:lar`, `:rfitzp`, `:toroidal` - `ising` -- singular-surface index for traceability + - `k_ref` -- reference-length ratio K = r_s·(dψ_N/dr) at the surface, + used by the Tearing runner to convert the ψ_N-referenced outer Δ' to + the r_s-referenced convention this layer works in (default `1.0`, + i.e. no conversion; `build_slayer_inputs` fills it from the equilibrium) + - `alpha_mercier` -- Mercier Frobenius exponent α = √(−D_I) for the same conversion + (default `0.5`, the slab/cylindrical value at D_I = −1/4) # Resistivity kwargs @@ -252,7 +271,9 @@ function slayer_parameters(; f_trap::Union{Real,Nothing}=nothing, nu_e_star::Union{Real,Nothing}=nothing, R_major_eff::Union{Real,Nothing}=nothing, - lnLambda_form::Symbol=:nrl) + lnLambda_form::Symbol=:nrl, + k_ref::Real=1.0, + alpha_mercier::Real=0.5) # Coulomb logarithm shared by the resistivity closure and τ_ee. lnLamb = coulomb_log_e(n_e, t_e; form=lnLambda_form) @@ -329,6 +350,16 @@ function slayer_parameters(; ) ) iota_e = Q_e / Q_e_minus_Q_i + iota_e > 0 || + throw( + ArgumentError( + "slayer_parameters: iota_e = Q_e/(Q_e - Q_i) = $iota_e is " * + "not positive, so the ion-sound-radius normalization " * + "D_norm = (d_beta/r_s)·S^(1/3)·sqrt(iota_e) is undefined. " * + "Check the diamagnetic-frequency inputs ω_*e, ω_*i — they " * + "are expected to carry opposite signs." + ) + ) # Plasma beta and compressibility lbeta = (5.0 / 3.0) * MU_0 * n_e * E_CHG * (t_e + t_i) / bt^2 @@ -342,7 +373,7 @@ function slayer_parameters(; # Normalized beta-related width and Δ-normalization d_beta = c_beta * d_i - D_norm = (d_beta / rs) * lu^(1.0 / 3.0) * sqrt(tau / (1.0 + tau)) + D_norm = (d_beta / rs) * lu^(1.0 / 3.0) * sqrt(iota_e) delta_n = lu^(1.0 / 3.0) / rs # Critical-Δ offset from chi_parallel matching @@ -360,6 +391,7 @@ function slayer_parameters(; rs=rs, R0=R0, bt=bt, sval_r=sval_r, dr_val=dr_val, dgeo_val=dgeo_val, eta=eta, d_beta=d_beta, - dc_tmp=dc_tmp, dc_type=dc_type + dc_tmp=dc_tmp, dc_type=dc_type, + k_ref=k_ref, alpha_mercier=alpha_mercier ) end diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index ece897c15..033f71eca 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -56,7 +56,7 @@ include("LayerInputs.jl") export SLAYERModel, SLAYERParameters, slayer_parameters export r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths -export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export surface_minor_radius, surface_da_dpsi, radial_label, build_slayer_inputs export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel end # module SLAYER diff --git a/src/KineticForces/CalculatedKineticMatrices.jl b/src/KineticForces/CalculatedKineticMatrices.jl index bb61135e4..d6436b710 100644 --- a/src/KineticForces/CalculatedKineticMatrices.jl +++ b/src/KineticForces/CalculatedKineticMatrices.jl @@ -37,6 +37,7 @@ section. - `ffs_intr`: ForceFreeStatesInternal (mode indexing) - `metric`: MetricData (provides ψ grid via `metric.xs`) - `mats`: MatrixSplines (used only for `numpert_total` cross-check) + - `psis`: ψ grid for the kernel; empty (default) evaluates on `metric.xs` # Keyword arguments @@ -59,9 +60,13 @@ function compute_calculated_kinetic_matrices( mats; kf_ctrl::KineticForcesControl=KineticForcesControl(), kinetic_profiles::Equilibrium.KineticProfileSplines, - species::Union{Nothing,AbstractVector{<:Equilibrium.ResolvedNTVSpecies}}=nothing + species::Union{Nothing,AbstractVector{<:Equilibrium.ResolvedNTVSpecies}}=nothing, + psis::Vector{Float64}=Float64[], + axis_psi_c::Float64=0.0 ) - xs = metric.xs + # The kernel is a pure function of psi (it evaluates equilibrium splines), so it can be + # driven over any knot list; default is the full equilibrium grid. + xs = isempty(psis) ? metric.xs : psis mpsi = length(xs) mpert = ffs_intr.mpert npert = ffs_intr.npert @@ -110,6 +115,20 @@ function compute_calculated_kinetic_matrices( splist = species === nothing ? [Equilibrium.ResolvedNTVSpecies(kf_ctrl.zi, kf_ctrl.mi, kf_ctrl.electron, "single", kinetic_profiles)] : species + + # Near-axis validity envelope: suppresses the drift-kinetic increments where the zero-orbit-width + # ordering fails (kernel evaluation skipped where it is 0). The same real scalar must scale kw and + # kt at the same ψ — that is what keeps the blocks Hermitian and the energy/torque relation intact. + env = ones(Float64, mpsi) + if kf_ctrl.axis_validity_suppression + psi_c = axis_psi_c + if psi_c > 0 + env .= kinetic_axis_validity_envelope.(xs, psi_c) + @info "Kinetic axis-validity suppression: psi_c=$(round(psi_c; sigdigits=3)), envelope reaches 1 at " * + "psi=$(round(2 * psi_c; sigdigits=3)) (kernel evaluation skipped below psi_c)" maxlog = 1 + end + end + for sp in splist # Hoisted per-species scalars: the closure then captures concrete values, not a # union-typed `sp` (the ternary above mixes the resolved vector with the fallback). @@ -124,6 +143,7 @@ function compute_calculated_kinetic_matrices( block_w = thread_block_w[tid] block_t = thread_block_t[tid] psi = xs[ipsi] + env[ipsi] == 0.0 && continue # inside the near-axis validity band: kernel skipped for in_idx in 1:npert n = ffs_intr.nlow + in_idx - 1 fill!(full_w, 0) @@ -148,8 +168,8 @@ function compute_calculated_kinetic_matrices( 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] + kw_flat[ipsi, idx, k] += env[ipsi] * full_w[i, j, k] + kt_flat[ipsi, idx, k] += env[ipsi] * full_t[i, j, k] end end end diff --git a/src/KineticForces/Compute.jl b/src/KineticForces/Compute.jl index 62c55d7f1..3ed4fdd9d 100644 --- a/src/KineticForces/Compute.jl +++ b/src/KineticForces/Compute.jl @@ -83,7 +83,7 @@ function integrate_psi_quadgk( wdfac::Float64, divxfac::Float64, electron::Bool, method::String, equil, intr::KineticForcesInternal, ctrl::KineticForcesControl, kinetic_profiles::Equilibrium.KineticProfileSplines; - psi_min::Float64=0.0, psi_max::Float64=1.0 + psi_min::Float64=0.0, psi_max::Float64=1.0, axis_psi_c::Float64=0.0 ) is_matrix_method = occursin("mm", method) mpert = intr.mpert @@ -100,6 +100,17 @@ function integrate_psi_quadgk( panel_psis=Float64[], resonance_psis=Float64[]) end + # Near-axis validity suppression on the torque density. `axis_psi_c` is the run's single + # boundary (widest-orbit species, cleared of rationals), so the quadrature domain matches the + # EL kinetic matrices and the Validity output; below it the integrand is identically zero. + psi_c = ctrl.axis_validity_suppression ? axis_psi_c : 0.0 + x0 = max(x0, psi_c) + 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[]) + end + psi_c > 0 && @info "Kinetic axis-validity suppression in ψ torque quadrature: domain starts at psi_c=$(round(psi_c; sigdigits=3))" + # The outer ψ-integral (the QuadGK batch / ψ-node loop) stays serial: QuadGK's refine # loop invokes the callback with small batches (~15 nodes), so threading it is # fork-join-bound. Instead thread the inner bounce-harmonic loop (2·nl+1 harmonics), @@ -149,6 +160,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm total += harm_vals[ell_idx] end + total *= kinetic_axis_validity_envelope(psi, psi_c) y[k] = total push!(logged_psi, psi) @@ -158,6 +170,7 @@ function integrate_psi_quadgk( for ell_idx in 1:nharm elems_accum .+= harm_elems[ell_idx] end + elems_accum .*= kinetic_axis_validity_envelope(psi, psi_c) push!(logged_elems, elems_accum) end end @@ -166,7 +179,7 @@ function integrate_psi_quadgk( # Panels at the rational surfaces the run resolved plus the kinetic-resonance # surfaces (thermal-energy Ω_ℓ = 0 for ℓ ∈ -nl:nl) — both are torque-density peaks. resonance_psis = kinetic_resonance_psi_nodes(kinetic_profiles, equil; n, nl, zi, mi, electron, wdfac) - pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis), x0, xout) + pts = psi_panel_points(vcat(intr.sing_psis, resonance_psis, [2 * psi_c]), x0, xout) bi = QuadGK.BatchIntegrand(psi_batch!, ComplexF64[], Float64[]) total, quad_err = quadgk(bi, pts...; atol=ctrl.atol_psi, rtol=ctrl.rtol_psi, maxevals=ctrl.maxevals_psi) @@ -291,10 +304,13 @@ block-diagonal kinetic matrices. - `ctrl::KineticForcesControl`: Control parameters specifying which methods to run - `equil`: PlasmaEquilibrium with 2D interpolants - `kinetic_profiles::Equilibrium.KineticProfileSplines`: Named kinetic-profile splines + - `axis_psi_c`: the run's near-axis validity boundary, suppressing the torque density below it + (`KineticForces.axis_validity_boundary`); 0 disables the suppression """ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticForcesInternal, ctrl::KineticForcesControl, equil, - kinetic_profiles::Equilibrium.KineticProfileSplines) + kinetic_profiles::Equilibrium.KineticProfileSplines; + axis_psi_c::Float64=0.0) for entry in METHOD_REGISTRY getfield(ctrl, entry.flag) || continue @@ -337,7 +353,7 @@ function compute_torque_all_methods!(state::KineticForcesState, intr::KineticFor n, ctrl.nl, ctrl.zi, ctrl.mi, ctrl.wdfac, ctrl.divxfac, ctrl.electron, method, equil, intr, ctrl, kinetic_profiles; - psi_min=ctrl.psilims[1], psi_max=ctrl.psilims[2]) + psi_min=ctrl.psilims[1], psi_max=ctrl.psilims[2], axis_psi_c=axis_psi_c) total_torque += result.total psi_nsteps_total += result.psi_nsteps diff --git a/src/KineticForces/KineticForcesStructs.jl b/src/KineticForces/KineticForcesStructs.jl index 0b58da3ec..16d895ada 100644 --- a/src/KineticForces/KineticForcesStructs.jl +++ b/src/KineticForces/KineticForcesStructs.jl @@ -3,12 +3,13 @@ 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. @@ -77,6 +78,7 @@ User-facing control parameters from the TOML `[KineticForces]` section. Configures which NTV methods to run, species parameters, tolerances, and output options. Constructed via keyword arguments or from a TOML dict: + ```julia ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForces"])...) ``` @@ -84,6 +86,13 @@ ctrl = KineticForcesControl(; (Symbol(k) => v for (k, v) in inputs["KineticForce Immutable: vary a field by building a new control rather than assigning to one (the multi-species loop does this per species, and `check_psi_quadrature_convergence`'s test builds a second control for its differing tolerance). + + - `axis_validity_suppression::Bool` - Suppress the calculated kinetic terms where the + zero-orbit-width ordering fails near the axis: below ψ_c — the outermost ψ at which a thermal + orbit width of the widest-orbit species reaches ⟨r⟩, moved clear of any rational surface the + transition would otherwise cut through — the increments are zeroed and the kernel is skipped, + rising to full strength at 2ψ_c through a C² envelope. Boundary and envelope come from the + equilibrium and kinetic profiles; there are no tuning parameters. Default `true`. """ @kwdef struct KineticForcesControl # Moment type @@ -121,32 +130,19 @@ builds a second control for its differing tolerance). nn::Int = 1 # Toroidal mode number nl::Int = 1 # Bounce harmonic number - # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy). - # Each level must be resolved more tightly than the one enclosing it, or the outer - # integrator chases its integrand's own quadrature noise instead of converging. - # *_xlmda: tolerances for the λ (pitch) integration - # *_x: tolerances for the x (energy) integration nested inside it; NaN ⇒ derive as - # nested_tolerance_margin × the pitch tolerances - # *_psi: tolerances for the outer ψ quadrature + # Tolerances, outermost to innermost: ψ quadrature ⊃ λ (pitch) ⊃ x (energy); each level must + # resolve tighter than the one enclosing it, or the outer integrator chases quadrature noise. atol_xlmda::Float64 = 1e-8 # Absolute tolerance for the inner pitch integration rtol_xlmda::Float64 = 1e-5 # Relative tolerance for the inner pitch integration atol_x::Float64 = NaN # Absolute tolerance for the energy integration (NaN ⇒ derived) rtol_x::Float64 = NaN # Relative tolerance for the energy integration (NaN ⇒ derived) - # The pitch integrand IS the energy integral, so the energy level is resolved this much - # tighter than the pitch level by default. + # Energy tolerances default to this fraction of the pitch ones (the pitch integrand is itself an energy integral). nested_tolerance_margin::Float64 = 1e-2 # Factor relating derived energy tolerances to the pitch ones - # rtol_psi is the primary convergence knob: ~2 significant figures matches the validity - # of the NTV model approximations. Do not set it tighter than the noise floor of the - # inner integrals (keep rtol_psi ≳ 10 × rtol_xlmda). + # Primary convergence knob (~2 sig figs matches NTV model validity); keep rtol_psi ≳ 10 × rtol_xlmda or it chases inner-integral noise. rtol_psi::Float64 = 1e-2 # Relative tolerance for outer ψ quadrature - # atol_psi is in N·m and therefore amplitude-sensitive: NTV scales as δB², so a 10× - # weaker applied field gives a 100× smaller torque and any fixed absolute tolerance can - # silently dominate termination with O(1) relative error. Default 0 (rtol-only); a - # nonzero value is an expert opt-out for near-net-zero-torque cases and triggers a - # warning when it dominates. + # Amplitude-sensitive (NTV ∝ δB²); default 0 (rtol-only). Nonzero is an expert opt-out for near-zero-torque cases and warns when it dominates. atol_psi::Float64 = 0.0 # Absolute tolerance for outer ψ quadrature [N·m] - # Runaway guard for the outer quadrature (e.g. sign-cancelling torque density with tiny - # net torque under rtol-only control); a convergence warning fires when hit. + # Runaway guard for rtol-only control (e.g. sign-cancelling torque density); warns when hit. maxevals_psi::Int = 2000 # Max integrand evaluations for outer ψ quadrature # Scaling factors @@ -158,6 +154,7 @@ builds a second control for its differing tolerance). nufac::Float64 = 1.0 # Collisionality scaling divxfac::Float64 = 1.0 # div(xi_perp) scaling + axis_validity_suppression::Bool = true # documented in the `## Fields` docstring # Energy integration parameters nutype::String = "harmonic" # Collision operator: "zero", "small", "krook", "harmonic" @@ -211,12 +208,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 @@ -306,17 +304,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 @@ -326,20 +324,22 @@ 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::PerturbedEquilibrium.PerturbedEquilibriumState, - ffs::ForceFreeStates.ForceFreeStatesResult, - equil::Equilibrium.PlasmaEquilibrium, - metric::ForceFreeStates.MetricData) + ffs::ForceFreeStates.ForceFreeStatesResult, + equil::Equilibrium.PlasmaEquilibrium, + metric::ForceFreeStates.MetricData) # Copy mode numbers from FFS kf_intr.mlow = ffs.mlow kf_intr.mhigh = ffs.mhigh @@ -408,9 +408,9 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur 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) @@ -430,8 +430,8 @@ function set_perturbation_data!(kf_intr::KineticForcesInternal, pe_state::Pertur 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 @@ -458,9 +458,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) @@ -530,8 +530,8 @@ Accumulated results from all KineticForces computations. Written to gpec.h5 under the "KineticForces" 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 75f527ae2..867fafff4 100644 --- a/src/KineticForces/Output.jl +++ b/src/KineticForces/Output.jl @@ -13,15 +13,16 @@ then write to gpec.h5 in a single pass. Write KineticForces results to the "KineticForces" group in gpec.h5. # Arguments -- `h5file::HDF5.File`: Open HDF5 file handle -- `state::KineticForcesState`: Accumulated computation results -- `dVdpsi_spline`: Optional dV/dψ_N profile interpolant; when given, dV/dψ_N is - written at the quadrature points so the torque density dT/dV = (dT/dψ)/(dV/dψ) - is directly available -- `species_label`: `nothing` writes the run total to `KineticForces//`; a label - (e.g. `"ion_z1_m2"`, `"electron"`) writes one species' contribution to - `KineticForces/PerSpecies/