diff --git a/docs/src/citations.md b/docs/src/citations.md index 014a258393..28a1292420 100644 --- a/docs/src/citations.md +++ b/docs/src/citations.md @@ -296,3 +296,19 @@ Provides the complete theory and implementation details, including NTV in the pr > DOI: [10.1063/1.4977898](https://doi.org/10.1063/1.4977898) Describes the fully self-consistent coupling between the perturbed equilibrium and neoclassical toroidal viscosity (NTV), providing the theoretical foundation for the KineticForces functionality. + +--- + +## Utilities + +> A. Redl, C. Angioni, E. Belli, O. Sauter, "A new set of analytical formulae for the computation of the bootstrap current and the neoclassical conductivity in tokamaks," +> *Physics of Plasmas* **28**, 022502 (2021). + +Source of the bootstrap-current coefficients L31, L32 and α (Eqs. 10-21) used by the `BootstrapCurrent` utilities, implemented as OpenFUSIONToolkit TokaMaker evaluates them. + +--- + +> O. Sauter, C. Angioni, and Y. R. Lin-Liu, "Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime," +> *Physics of Plasmas* **6**, 2834 (1999). + +Defines the circulating and trapped particle fractions and the collisionality conventions used by `BootstrapCurrent` and `NeoclassicalResistivity`. diff --git a/docs/src/utilities.md b/docs/src/utilities.md index 3eaf4fa987..8233cfd515 100644 --- a/docs/src/utilities.md +++ b/docs/src/utilities.md @@ -23,6 +23,17 @@ neoclassical models) used to set the Lundquist number in the tearing stack. Modules = [GeneralizedPerturbedEquilibrium.Utilities.NeoclassicalResistivity] ``` +## Bootstrap current + +The flux-surface-averaged bootstrap current ``\langle j_{bs}\cdot B\rangle`` from the Redl et al. (2021) closure, assembled as OpenFUSIONToolkit TokaMaker assembles it: exact Redl coefficients with ``L_{34}=L_{31}``, the Koh ion collisionality, NRL Coulomb logarithms floored at 10, and the Sauter circulating fraction from surface moments of ``|B|``. +ψ is the poloidal flux per radian, increasing from the axis to the edge, so ``\psi = \psi_o\hat\psi``; ``I(\psi)=RB_\phi`` is signed and the result is co-directed with the equilibrium ``\langle j\cdot B\rangle`` of the same convention. +Temperatures are in eV and densities in m⁻³. +See Redl et al., Phys. Plasmas **28**, 022502 (2021) and, for the circulating fraction, Sauter et al., Phys. Plasmas **6**, 2834 (1999). + +```@autodocs +Modules = [GeneralizedPerturbedEquilibrium.Utilities.BootstrapCurrent] +``` + ## HDF5 annotations Self-describing metadata for `gpec.h5` (long_name/units/dims attributes and HDF5 diff --git a/src/Utilities/BootstrapCurrent.jl b/src/Utilities/BootstrapCurrent.jl new file mode 100644 index 0000000000..ea015e0726 --- /dev/null +++ b/src/Utilities/BootstrapCurrent.jl @@ -0,0 +1,235 @@ +""" + BootstrapCurrent + +Flux-surface-averaged bootstrap current ⟨j_bs·B⟩ from the Redl et al. 2021 closure, assembled the way +OpenFUSIONToolkit TokaMaker assembles it: exact Redl coefficients with L34 = L31, Koh ion collisionality, +NRL electron and ion Coulomb logarithms floored at 10, and the Sauter circulating fraction from surface moments of |B|. + +Conventions: ψ is the poloidal flux per radian [Wb/rad] increasing from the axis to the edge, so ψ = psio ψ_N and the +gradients passed to [`jdotB_bootstrap`](@ref) are d/dψ of profiles that fall outward. I(ψ) = R B_φ is signed and +⟨j_bs·B⟩ = −I p_e [L31 p′/p_e + L32 T_e′/T_e + L34 α (1 − R_pe)/R_pe T_i′/T_i] is co-directed with the equilibrium +⟨j·B⟩ = −(F p′ + F′⟨B²⟩/μ₀) of the same convention. Temperatures are in eV, densities in m⁻³, lengths in m. + +References: Redl et al., Phys. Plasmas 28, 022502 (2021); Sauter et al., Phys. Plasmas 6, 2834 (1999). +""" +module BootstrapCurrent + +using FastInterpolations +using ..NeoclassicalResistivity: coulomb_log_e + +export redl_coefficients, coulomb_log_i, nu_star_i_koh, nu_star_e_floored, circulating_fraction, jdotB_bootstrap, surface_moments, bootstrap_current + +# Elementary charge as used by TokaMaker (CODATA 2018) +const E_CHARGE = 1.602176634e-19 + +""" + redl_coefficients(f_t, Z_eff, nu_e_star, nu_i_star) -> (; L31, L32, L34, alpha) + +Redl 2021 Eqs. 10-21 with L34 = L31, as in TokaMaker `redl_bootstrap`. +""" +function redl_coefficients(f_t::Real, Z_eff::Real, nu_e_star::Real, nu_i_star::Real) + fT, Z, nue, nui = f_t, Z_eff, nu_e_star, nu_i_star + sZ1 = sqrt(max(Z - 1, 0.0)) + # L31 (Eqs. 10-11) + d1 = 0.67 * (1 - 0.7fT) * sqrt(nue) / (0.56 + 0.44Z) + d2 = (0.52 + 0.086 * sqrt(nue)) * (1 + 0.87fT) * nue / (1 + 1.13sZ1) + X31 = fT / (1 + d1 + d2) + dZ = Z^1.2 - 0.71 + L31 = (1 + 0.15 / dZ) * X31 - (0.22 / dZ) * X31^2 + (0.01 / dZ) * X31^3 + (0.06 / dZ) * X31^4 + # L32 (Eqs. 12-16) + dee2 = 0.23 * (1 - 0.96fT) * sqrt(nue) / sqrt(Z) + dee3 = 0.13 * (1 - 0.38fT) * nue / Z^2 * (sqrt(1 + 2sZ1) + fT^2 * sqrt((0.075 + 0.25 * (Z - 1)^2) * nue)) + Xee = fT / (1 + dee2 + dee3) + Fee = (0.1 + 0.6Z) / (Z * (0.77 + 0.63 * (1 + max(Z - 1, 0.0)^1.1))) * (Xee - Xee^4) + + 0.7 / (1 + 0.2Z) * (Xee^2 - Xee^4 - 1.2 * (Xee^3 - Xee^4)) + 1.3 / (1 + 0.5Z) * Xee^4 + dei2 = 0.87 * (1 + 0.39fT) * sqrt(nue) / (1 + 2.95 * (Z - 1)^2) + dei3 = 1.53 * (1 - 0.37fT) * nue * (2 + 0.375 * (Z - 1)) + Xei = fT / (1 + dei2 + dei3) + Fei = -(0.4 + 1.93Z) / (Z * (0.8 + 0.6Z)) * (Xei - Xei^4) + 5.5 / (1.5 + 2Z) * (Xei^2 - Xei^4 - 0.8 * (Xei^3 - Xei^4)) - 1.3 / (1 + 0.5Z) * Xei^4 + L32 = Fee + Fei + # α (Eqs. 20-21) + alpha0 = -(0.62 + 0.055 * (Z - 1)) / (0.53 + 0.17 * (Z - 1)) * (1 - fT) / (1 - (0.31 - 0.065 * (Z - 1)) * fT - 0.25fT^2) + alpha = ((alpha0 + 0.7Z * sqrt(fT) * sqrt(nui)) / (1 + 0.18 * sqrt(nui)) - 0.002 * nui^2 * fT^6) / (1 + 0.004 * nui^2 * fT^6) + return (; L31, L32, L34=L31, alpha) +end + +""" + coulomb_log_i(n_e, n_i, T_i) -> Float64 + +Ion-ion Coulomb logarithm `30 − ln(Z_avg³ √n_i / T_i^1.5)` with `Z_avg = max(n_e/n_i, 1)`, floored at 10 (TokaMaker form). +`n` in m⁻³, `T_i` in eV. +""" +function coulomb_log_i(n_e::Real, n_i::Real, T_i::Real) + Z_avg = max(n_e / n_i, 1.0) + return max(30.0 - log(Z_avg^3 * sqrt(n_i) / T_i^1.5), 10.0) +end + +""" + nu_star_e_floored(n_e, T_e, R_geo, eps, q, Z_eff) -> Float64 + +Electron collisionality of `nu_star_e` with the NRL Coulomb logarithm floored at 10, as TokaMaker. +""" +function nu_star_e_floored(n_e::Real, T_e::Real, R_geo::Real, eps::Real, q::Real, Z_eff::Real) + lnL = max(coulomb_log_e(n_e, T_e), 10.0) + return 6.921e-18 * abs(q) * R_geo * n_e * Z_eff * lnL / (T_e^2 * eps^1.5) +end + +""" + nu_star_i_koh(n_e, n_i, T_i, R_geo, eps, q, Z_eff; Z_dom=1) -> Float64 + +Ion collisionality, Koh form: `4.90e-18 |q| R n_i Z_ion⁴ lnΛ_ii / (T_i² ε^1.5)` with +`Z_ion = (Z_dom² Z_avg Z_eff)^(1/4)` and `Z_avg = n_e/n_i`. +""" +function nu_star_i_koh(n_e::Real, n_i::Real, T_i::Real, R_geo::Real, eps::Real, q::Real, Z_eff::Real; Z_dom::Real=1) + Z_ion4 = Z_dom^2 * (n_e / n_i) * Z_eff + return 4.90e-18 * abs(q) * R_geo * n_i * Z_ion4 * coulomb_log_i(n_e, n_i, T_i) / (T_i^2 * eps^1.5) +end + +""" + circulating_fraction(B, weights) -> Float64 + +Circulating fraction `f_c = 1 − (3/4 f_tu + 1/4 f_tl)` of Sauter 1999 as TokaMaker evaluates it from moments of |B| on a surface, +with `f_tu = 1 − ⟨B²⟩/⟨B⟩² (1 − √(1−h)(1+h/2))`, `f_tl = 1 − ⟨B²⟩/B_max² ⟨(1 − √(1−b)(1+b/2))/b²⟩`, `h = ⟨B⟩/B_max`, `b = B/B_max`. +`B` samples |B| on a periodic uniform θ grid (no repeated endpoint) and `weights` are the flux-surface-average weights (the Jacobian). +""" +function circulating_fraction(B::AbstractVector{<:Real}, weights::AbstractVector{<:Real}) + wsum = sum(weights) + Bmax = maximum(B) + avgB = sum(weights .* B) / wsum + avgB2 = sum(weights .* B .^ 2) / wsum + h = avgB / Bmax + h2 = avgB2 / Bmax^2 + ftu = 1 - h2 / h^2 * (1 - sqrt(1 - h) * (1 + h / 2)) + g(b) = (1 - sqrt(1 - b) * (1 + b / 2)) / b^2 + hf = sum(weights .* g.(min.(B ./ Bmax, 1.0))) / wsum + ftl = 1 - h2 * hf + return 1 - (0.75ftu + 0.25ftl) +end + +""" + jdotB_bootstrap(; T_e, T_i, n_e, n_i, Z_eff, f_t, nu_e_star, nu_i_star, I_psi, dT_e, dT_i, dn_e, dn_i) -> NamedTuple + +TokaMaker `jboot1` assembly of ⟨j_bs·B⟩ [T A/m²]. Temperatures in eV, densities in m⁻³, `I_psi = R B_φ` [T m], and the four +gradients are d/dψ with ψ in Wb/rad increasing outward. Returns `(; jdotB, L31, L32, L34, alpha)`. +""" +function jdotB_bootstrap(; T_e, T_i, n_e, n_i, Z_eff, f_t, nu_e_star, nu_i_star, I_psi, dT_e, dT_i, dn_e, dn_i) + c = redl_coefficients(f_t, Z_eff, nu_e_star, nu_i_star) + p_e = E_CHARGE * n_e * T_e + p_i = E_CHARGE * n_i * T_i + R_pe = p_e / (p_e + p_i) + dp = (n_e * dT_e + T_e * dn_e + n_i * dT_i + T_i * dn_i) * E_CHARGE + bra = c.L31 * dp / p_e + c.L32 * dT_e / T_e + c.L34 * c.alpha * (1 - R_pe) / R_pe * dT_i / T_i + return (; jdotB=-I_psi * p_e * bra, c.L31, c.L32, c.L34, c.alpha) +end + +""" + surface_moments(equil, psi_N) -> NamedTuple + +Geometry of the flux surface `psi_N` of a `PlasmaEquilibrium`, sampled on its θ grid: `B` (|B| per θ), `weights` (Jacobian), `R`, +and `q`, `F = R B_φ`, `R_max`, `R_min`, `R_geo = (R_max + R_min)/2`, `eps = (R_max − R_min)/(R_max + R_min)`, `avg_B2 = ⟨B²⟩`, +`avg_invR = ⟨1/R⟩`, `avg_invR2 = ⟨1/R²⟩`. Hamada-style metric as in the equilibrium construction; `profiles.F_spline` is 2πF. +""" +function surface_moments(equil, psi_N::Real) + profiles = equil.profiles + chi1 = 2π * equil.psio + psi = Float64(psi_N) + q = profiles.q_spline(psi) + F = profiles.F_spline(psi) / 2π + ys = equil.rzphi_ys[1:(end-1)] + n = length(ys) + B = zeros(n) + w = zeros(n) + R = zeros(n) + for (i, theta) in enumerate(ys) + rfac = sqrt(equil.rzphi_rsquared((psi, theta))) + jac = equil.rzphi_jac((psi, theta)) + fy1 = FastInterpolations.deriv_view(equil.rzphi_rsquared, (0, 1))((psi, theta)) + fy2 = FastInterpolations.deriv_view(equil.rzphi_offset, (0, 1))((psi, theta)) + fy3 = FastInterpolations.deriv_view(equil.rzphi_nu, (0, 1))((psi, theta)) + eta = 2π * (theta + equil.rzphi_offset((psi, theta))) + r = equil.ro + rfac * cos(eta) + v21 = fy1 / (2 * rfac * jac) + v22 = (1 + fy2) * 2π * rfac / jac + v23 = fy3 * r / jac + v33 = 2π * r / jac + B[i] = chi1 * sqrt(v21^2 + v22^2 + (v23 + q * v33)^2) + w[i] = jac + R[i] = r + end + wsum = sum(w) + avg(x) = sum(w .* x) / wsum + R_max, R_min = maximum(R), minimum(R) + return (; B, weights=w, R, q, F, R_max, R_min, R_geo=(R_max + R_min) / 2, eps=(R_max - R_min) / (R_max + R_min), + avg_B2=avg(B .^ 2), avg_invR=avg(1 ./ R), avg_invR2=avg(1 ./ R .^ 2)) +end + +# Second-order differences on a non-uniform grid, as TokaMaker's gradient_ (numpy.gradient, edge_order=2) +function _gradient(x::AbstractVector, y::AbstractVector) + n = length(x) + d = zeros(n) + for i in 2:(n-1) + hl, hr = x[i] - x[i-1], x[i+1] - x[i] + d[i] = (hl^2 * y[i+1] - hr^2 * y[i-1] + (hr^2 - hl^2) * y[i]) / (hl * hr * (hl + hr)) + end + h1, h2 = x[2] - x[1], x[3] - x[2] + d[1] = -(2h1 + h2) / (h1 * (h1 + h2)) * y[1] + (h1 + h2) / (h1 * h2) * y[2] - h1 / (h2 * (h1 + h2)) * y[3] + h1, h2 = x[n-1] - x[n-2], x[n] - x[n-1] + d[n] = h2 / (h1 * (h1 + h2)) * y[n-2] - (h1 + h2) / (h1 * h2) * y[n-1] + (2h2 + h1) / (h2 * (h1 + h2)) * y[n] + return d +end + +""" + bootstrap_current(equil, kinetic=equil.kinetic; psi_N=range(0.05, 0.95; length=91), Z_dom=1) -> NamedTuple + +Bootstrap current profiles on `psi_N` for a `PlasmaEquilibrium` and its `KineticProfileData` (needs `n_e`, `n_i`, `T_e`, `T_i`; `Z_eff` is the +`zeff` profile, else quasi-neutrality with one impurity of `impurity_charge` and a charge-1 main ion). The kinetic columns are cubic-spline interpolated onto `psi_N` and differentiated by +second-order finite differences in ψ = psio ψ_N, as TokaMaker does, so a fine `psi_N` grid is needed for smooth gradients. +`psi_N` must hold at least three strictly increasing points inside (0, 1). +Returns `psi_N`, `jdotB` (⟨j_bs·B⟩ [T A/m²]), `jphi_bs = ⟨j_bs·B⟩ F⟨1/R⟩/⟨B²⟩` [A/m²], `f_t`, `nu_e_star`, `nu_i_star`, `eps`, `R_geo`, `q`, `F`, +`avg_B2`, `avg_invR`, `avg_invR2`, `L31`, `L32`, `L34`, `alpha`. `jphi_bs` is the field-aligned part only; TokaMaker's `jphi_bs` +also adds the pressure-driven term P′(⟨R⟩ − F²⟨1/R⟩/⟨B²⟩). +""" +function bootstrap_current(equil, kinetic=equil.kinetic; psi_N=range(0.05, 0.95; length=91), Z_dom::Real=1) + kinetic === nothing && throw(ArgumentError("bootstrap_current: no kinetic profiles attached")) + for name in (:n_e, :n_i, :T_e, :T_i) + getfield(kinetic, name) === nothing && throw(ArgumentError("bootstrap_current: kinetic profile $name is missing")) + end + x = collect(Float64, psi_N) + (length(x) >= 3 && all(>(0), diff(x)) && 0 < x[1] && x[end] < 1) || + throw(ArgumentError("bootstrap_current: psi_N needs at least three strictly increasing points inside (0, 1)")) + interp(y) = (itp = cubic_interp(kinetic.psi, y); [itp(xi) for xi in x]) + ne, ni, Te, Ti = interp(kinetic.n_e), interp(kinetic.n_i), interp(kinetic.T_e), interp(kinetic.T_i) + # Z_eff as the inner layers take it: the file's profile, else quasi-neutrality with one impurity and a charge-1 main ion; + # the floor at 1 hides an n_i > n_e inconsistency in the file. + # An all-zero `zeff` column counts as absent. + Zeff = if kinetic.zeff !== nothing && !all(iszero, kinetic.zeff) + minimum(kinetic.zeff) >= 1 - 1e-8 || throw(ArgumentError("bootstrap_current: zeff = $(minimum(kinetic.zeff)) < 1 is unphysical")) + interp(kinetic.zeff) + else + zimp = Float64(kinetic.impurity_charge) + zimp >= 1 || throw(ArgumentError("bootstrap_current: impurity_charge = $zimp < 1 is unphysical")) + [max(zimp - (nii / nee) * (zimp - 1), 1.0) for (nii, nee) in zip(ni, ne)] + end + psi = equil.psio .* x + dne, dni, dTe, dTi = _gradient(psi, ne), _gradient(psi, ni), _gradient(psi, Te), _gradient(psi, Ti) + geo = [surface_moments(equil, xi) for xi in x] + n = length(x) + out = (; psi_N=x, jdotB=zeros(n), jphi_bs=zeros(n), f_t=zeros(n), nu_e_star=zeros(n), nu_i_star=zeros(n), eps=[g.eps for g in geo], R_geo=[g.R_geo for g in geo], + q=[g.q for g in geo], F=[g.F for g in geo], avg_B2=[g.avg_B2 for g in geo], avg_invR=[g.avg_invR for g in geo], avg_invR2=[g.avg_invR2 for g in geo], + L31=zeros(n), L32=zeros(n), L34=zeros(n), alpha=zeros(n)) + for i in 1:n + g = geo[i] + f_t = 1 - circulating_fraction(g.B, g.weights) + nue = nu_star_e_floored(ne[i], Te[i], g.R_geo, g.eps, g.q, Zeff[i]) + nui = nu_star_i_koh(ne[i], ni[i], Ti[i], g.R_geo, g.eps, g.q, Zeff[i]; Z_dom) + r = jdotB_bootstrap(; T_e=Te[i], T_i=Ti[i], n_e=ne[i], n_i=ni[i], Z_eff=Zeff[i], f_t, nu_e_star=nue, nu_i_star=nui, I_psi=g.F, + dT_e=dTe[i], dT_i=dTi[i], dn_e=dne[i], dn_i=dni[i]) + out.jdotB[i] = r.jdotB + out.jphi_bs[i] = r.jdotB * g.F * g.avg_invR / g.avg_B2 + out.f_t[i], out.nu_e_star[i], out.nu_i_star[i] = f_t, nue, nui + out.L31[i], out.L32[i], out.L34[i], out.alpha[i] = r.L31, r.L32, r.L34, r.alpha + end + return out +end + +end # module BootstrapCurrent diff --git a/src/Utilities/Utilities.jl b/src/Utilities/Utilities.jl index 466783d3fb..1651ec8467 100644 --- a/src/Utilities/Utilities.jl +++ b/src/Utilities/Utilities.jl @@ -13,6 +13,7 @@ mathematical utilities. - `PhysicalConstants`: SI physical constants matching Fortran GPEC/SLAYER values - `NeoclassicalResistivity`: Spitzer/Sauter/Redl resistivity closures shared by the GGJ and SLAYER inner-layer models + - `BootstrapCurrent`: Redl bootstrap current ⟨j_bs·B⟩ as OpenFUSIONToolkit TokaMaker computes it - `HDF5Annotations`: self-describing metadata (long_name/units/dims attributes and HDF5 Dimension Scales) for the gpec.h5 output - `HermiteFEM`: the Hermite-cubic basis and packed grid shared by the Galerkin solvers @@ -24,6 +25,7 @@ include("FourierCoefficients.jl") include("PhysicalConstants.jl") include("KineticProfiles.jl") include("NeoclassicalResistivity.jl") +include("BootstrapCurrent.jl") include("HDF5Annotations.jl") include("Output.jl") include("HDF5Schema.jl") @@ -56,4 +58,8 @@ export coulomb_log_e, eta_spitzer, tau_ee_spitzer_harm, eta_spitzer_harm export trapped_fraction, trapped_fraction_eps export nu_star_e, eta_neoclassical +using .BootstrapCurrent +export BootstrapCurrent +export redl_coefficients, circulating_fraction, jdotB_bootstrap, bootstrap_current + end # module Utilities diff --git a/test/runtests_bootstrap_current.jl b/test/runtests_bootstrap_current.jl new file mode 100644 index 0000000000..956ff311eb --- /dev/null +++ b/test/runtests_bootstrap_current.jl @@ -0,0 +1,121 @@ +using Test +using TOML +using GeneralizedPerturbedEquilibrium + +@testset "Bootstrap current" begin + GPEC = GeneralizedPerturbedEquilibrium + BC = GPEC.Utilities.BootstrapCurrent + Eq = GPEC.Equilibrium + EC = 1.602176634e-19 + + # (f_t, Z_eff, nu_e*, nu_i*, L31, L32, alpha, jdotB) from OpenFUSIONToolkit TokaMaker `redl_bootstrap` (commit 637ddb5), jboot1, use_legacy_L34=False + redl_oracle = [ + (0.1, 1.0, 0.01, 0.01, 0.13545033708131518, -0.08305470624423099, -1.0483168453295013, 53865.33136495606), + (0.25, 1.5, 0.1, 0.05, 0.2359363283734479, -0.0715742719918062, -0.7167460948754181, 130260.29898267234), + (0.4, 2.0, 1.0, 0.5, 0.24743537097926352, 0.06276713461942436, -0.04000751746519661, 209222.55685617123), + (0.55, 1.2, 3.0, 1.0, 0.184850306311777, 0.14119356432215222, -0.03275606630693376, 183761.98078772676), + (0.7, 2.5, 10.0, 5.0, 0.10564861564801808, 0.151141790537118, 2.0418011480760123, 168358.79959828986), + (0.35, 3.0, 0.5, 10.0, 0.27010731762521395, 0.020969260841897194, 2.1182499511424036, 329096.2052590009) + ] + # (n_e, n_i, T_e, T_i, lnLambda_e, lnLambda_ii) from TokaMaker `calculate_ln_lambda` (NRL electron, Zavg ion) + lnl_oracle = [ + (4e+19, 3.5e+19, 1500.0, 1200.0, 15.653266337878941, 17.73358120809945), + (1e+20, 9e+19, 5000.0, 4000.0, 16.399094471655047, 19.151822241068018), + (2e+19, 2e+19, 300.0, 300.0, 14.390400380045524, 16.334541738260896), + (5e+18, 4e+18, 20000.0, 18000.0, 19.283255507827988, 22.611346883868467) + ] + + @testset "Redl coefficients and assembly match TokaMaker" begin + for (ft, Z, nue, nui, L31, L32, alpha, jB) in redl_oracle + c = BC.redl_coefficients(ft, Z, nue, nui) + @test c.L31 ≈ L31 rtol = 1e-12 + @test c.L32 ≈ L32 rtol = 1e-12 + @test c.L34 ≈ L31 rtol = 1e-12 + @test c.alpha ≈ alpha rtol = 1e-12 + Te, Ti, ne, ni = 1500.0, 1200.0, 4.0e19, 3.5e19 + r = BC.jdotB_bootstrap(; T_e=Te, T_i=Ti, n_e=ne, n_i=ni, Z_eff=Z, f_t=ft, nu_e_star=nue, nu_i_star=nui, I_psi=5.0, + dT_e=-9000.0, dT_i=-7000.0, dn_e=-1.5e20, dn_i=-1.2e20) + @test r.jdotB ≈ jB rtol = 1e-12 + end + end + + @testset "Coulomb logarithms and collisionality" begin + for (ne, ni, Te, Ti, le, li) in lnl_oracle + @test max(BC.coulomb_log_e(ne, Te), 10.0) ≈ le rtol = 1e-12 + @test BC.coulomb_log_i(ne, ni, Ti) ≈ li rtol = 1e-12 + end + # Floors at 10 + @test BC.coulomb_log_i(1e26, 1e26, 1.0) == 10.0 + # Koh form with Z_ion⁴ = Z_dom² (n_e/n_i) Z_eff, written out + ne, ni, Te, Ti, R, eps, q, Z = 4e19, 3.5e19, 1500.0, 1200.0, 1.7, 0.3, 2.0, 1.8 + lnii = 30 - log((ne / ni)^3 * sqrt(ni) / Ti^1.5) + @test BC.nu_star_i_koh(ne, ni, Ti, R, eps, q, Z) ≈ 4.90e-18 * q * R * ni * (ne / ni) * Z * lnii / (Ti^2 * eps^1.5) rtol = 1e-12 + @test BC.nu_star_i_koh(ne, ni, Ti, R, eps, -q, Z; Z_dom=2) ≈ 4 * BC.nu_star_i_koh(ne, ni, Ti, R, eps, q, Z) rtol = 1e-12 + lne = max(23.5 - log(sqrt(ne / 1e6) * Te^-1.25) - sqrt(1e-5 + (log(Te) - 2)^2 / 16), 10.0) + @test BC.nu_star_e_floored(ne, Te, R, eps, q, Z) ≈ 6.921e-18 * q * R * ne * Z * lne / (Te^2 * eps^1.5) rtol = 1e-12 + end + + @testset "Circulating fraction" begin + # Model surface B = B0/(1 + ε cos θ) with uniform weights; the reference evaluates the definition term by term with plain loops + function reference_fc(eps, N) + thetas = [2π * (k - 1) / N for k in 1:N] + Bs = [1 / (1 + eps * cos(t)) for t in thetas] + Bmax = maximum(Bs) + h = sum(Bs) / N / Bmax + h2 = sum(Bs .^ 2) / N / Bmax^2 + ftu = 1 - h2 / h^2 * (1 - sqrt(1 - h) * (1 + h / 2)) + hf = sum((1 - sqrt(1 - Bi / Bmax) * (1 + Bi / Bmax / 2)) / (Bi / Bmax)^2 for Bi in Bs) / N + return 1 - (0.75ftu + 0.25 * (1 - h2 * hf)) + end + n = 512 + th = range(0, 2π; length=n + 1)[1:n] + for eps in (0.05, 0.2, 0.5) + B = 1 ./ (1 .+ eps .* cos.(th)) + @test BC.circulating_fraction(B, ones(n)) ≈ reference_fc(eps, n) rtol = 1e-12 + # Grid convergence against a 100x finer grid; the √(1−b) cusp at B_max limits the periodic trapezoid to O(h²) + @test BC.circulating_fraction(B, ones(n)) ≈ reference_fc(eps, 100n) rtol = 1e-5 + end + # Small-ε limit f_t → 1.46 √ε + eps = 1e-4 + B = 1 ./ (1 .+ eps .* cos.(th)) + @test 1 - BC.circulating_fraction(B, ones(n)) ≈ 1.46 * sqrt(eps) rtol = 0.1 + # The 3/4-1/4 Sauter form approximates the exact f_c = 3/4 ⟨B²⟩ ∫_0^{1/Bmax} λ dλ / ⟨√(1 − λB)⟩ to about a percent + for eps in (0.1, 0.3) + B = 1 ./ (1 .+ eps .* cos.(th)) + lam = ((1:20000) .- 0.5) ./ 20000 ./ maximum(B) + integral = sum(l / (sum(sqrt.(1 .- l .* B)) / n) for l in lam) * step(lam) + fc_exact = 0.75 * sum(B .^ 2) / n * integral + @test BC.circulating_fraction(B, ones(n)) ≈ fc_exact rtol = 0.02 + end + end + + @testset "Second-order gradient is exact for quadratics on a non-uniform grid" begin + x = [0.0, 0.1, 0.25, 0.3, 0.5, 0.8, 0.85, 1.0] + @test BC._gradient(x, 3 .* x .^ 2 .- 2 .* x .+ 1) ≈ 6 .* x .- 2 rtol = 1e-12 + end + + @testset "End to end on a DIII-D-like equilibrium with kinetic profiles" begin + dir = joinpath(@__DIR__, "..", "examples", "DIIID-like_SLAYER_example") + inputs = TOML.parsefile(joinpath(dir, "gpec.toml")) + equil = Eq.PlasmaEquilibrium(GPEC.equilibrium_config(inputs, dir)) + x = collect(range(0.05, 0.95; length=91)) + r = BC.bootstrap_current(equil; psi_N=x) + @test all(isfinite, r.jdotB) + @test all(0 .< r.f_t .< 1) + @test all(r.eps .> 0) + # Equilibrium ⟨j·B⟩ = −(F p′ + F′⟨B²⟩/μ0) and toroidal ⟨j_φ/R⟩ = −(p′ + F F′⟨R⁻²⟩/μ0) in the same ψ convention (profiles hold μ0 dp/dψ_N and 2πF) + prof, psio = equil.profiles, equil.psio + jdotB_eq = [-(r.F[i] * prof.P_deriv(xi) / psio + prof.F_deriv(xi) / (2π) / psio * r.avg_B2[i]) for (i, xi) in enumerate(x)] ./ (4π * 1e-7) + @test sign(sum(r.jdotB)) == sign(sum(jdotB_eq)) + # Bootstrap fraction of the toroidal current, I = ∫⟨j_φ/R⟩ dV/2π + jphiR_bs = r.jdotB .* r.F .* r.avg_invR2 ./ r.avg_B2 + jphiR_tot = [-(prof.P_deriv(xi) + r.F[i] * prof.F_deriv(xi) / (2π) * r.avg_invR2[i]) / psio for (i, xi) in enumerate(x)] ./ (4π * 1e-7) + vp = [prof.dVdpsi_spline(xi) for xi in x] + trapz(f) = sum((f[1:(end-1)] .+ f[2:end]) ./ 2 .* diff(x)) + frac = trapz(jphiR_bs .* vp) / trapz(jphiR_tot .* vp) + @test 0.02 < frac < 0.8 + @test_throws ArgumentError BC.bootstrap_current(equil; psi_N=[0.2, 0.5]) + @test_throws ArgumentError BC.bootstrap_current(equil; psi_N=[0.5, 0.3, 0.7]) + @test_throws ArgumentError BC.bootstrap_current(equil; psi_N=[0.0, 0.5, 0.9]) + end +end diff --git a/test/suites.jl b/test/suites.jl index dfe79f9e4f..cbf8601b41 100644 --- a/test/suites.jl +++ b/test/suites.jl @@ -1,6 +1,7 @@ # Keep the full suite order while assigning each file to exactly one CI group. const TEST_SUITES = [ "runtests_utilities.jl" => "core", + "runtests_bootstrap_current.jl" => "core", "runtests_fouriertransforms.jl" => "core", "runtests_vacuum.jl" => "core", "runtests_equil.jl" => "core",