diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index bee2a8a5e..8c8c31012 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -24,6 +24,7 @@ import .GGJ: solve_inner_converged # experimental, not exported (reachable as a 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: set_resistive_width_based_psihigh, print_layer_width_table export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner export GGJ, GGJModel, GGJParameters @@ -34,5 +35,6 @@ export glasser_wang_2020_eq55 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 set_resistive_width_based_psihigh, print_layer_width_table end # module InnerLayer diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl index eb16a7639..d59f5273e 100644 --- a/src/InnerLayer/SLAYER/LayerThickness.jl +++ b/src/InnerLayer/SLAYER/LayerThickness.jl @@ -137,6 +137,8 @@ is built from, retained as a drift-scale reference. - `delta_s` -- complex layer thickness `δ_s = dels_db · d_β` [m] - `delta_s_m` -- `|δ_s|`, the resistive layer thickness [m] (primary) - `d_beta` -- β-weighted ion scale `c_β·d_i` [m] (drift reference) + - `δ_FKR` -- FKR resistive-inertial length scale, `r_s · S^(-1/3)` [m] + - `δ_visco` -- visco-resistive length scale, `δ_FKR · P_perp^(1/6)` [m] `delta_s_m` should sit within a few orders of magnitude of `d_beta` for a well-posed surface (`dels_db` is O(1)); a large gap flags a normalisation @@ -150,6 +152,9 @@ struct LayerWidths delta_s::ComplexF64 delta_s_m::Float64 d_beta::Float64 + δ_FKR::Float64 + δ_visco::Float64 + end """ @@ -165,7 +170,9 @@ to `riccati_del_s`. function slayer_layer_thickness(p::SLAYERParameters; kwargs...) dels_db = riccati_del_s(p; kwargs...) delta_s = dels_db * p.d_beta + δ_FKR = p.rs * p.lu^(-1.0/3.0) + δ_visco = δ_FKR * p.P_perp^(1.0/6.0) return LayerWidths(p.ising, p.m, p.n, dels_db, delta_s, abs(delta_s), - p.d_beta) + p.d_beta, δ_FKR, δ_visco) end diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index ece897c15..0e9a80dcd 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -21,6 +21,8 @@ module SLAYER using LinearAlgebra using StaticArrays +using FastInterpolations +using HDF5 import ..InnerLayerModel, ..InnerLayerResponse, ..solve_inner, ..InnerLayerParameters using ...Utilities.PhysicalConstants @@ -52,11 +54,13 @@ include("LayerParameters.jl") include("Riccati.jl") include("LayerThickness.jl") include("LayerInputs.jl") +include("SetResistiveWidthPsihigh.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 set_resistive_width_based_psihigh, print_layer_width_table export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel end # module SLAYER diff --git a/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl b/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl new file mode 100644 index 000000000..5ca01ae71 --- /dev/null +++ b/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl @@ -0,0 +1,158 @@ +function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0, z_i=1.0) + equil = r.equil + + x_knots = equil.profiles.xs + q_knots = equil.profiles.q_spline.y + iota_itp_real = cubic_interp(x_knots, 1.0 ./ q_knots; extrap=Extrap(:extend)) + q_of_psi_real = ψ -> 1.0 / iota_itp_real(ψ) + dqdpsi_real = ψ -> -iota_itp_real(ψ; deriv=DerivOp(1)) / iota_itp_real(ψ)^2 + + ψ1, ψ2 = x_knots[end-1], x_knots[end] + q1, q2 = q_of_psi_real(ψ1), q_of_psi_real(ψ2) + A = (q2 - q1) / (log(1-ψ1) - log(1-ψ2)) + ψ_last, q_last = ψ2, q2 + + q_of_psi_hybrid(ψ) = ψ <= ψ_last ? q_of_psi_real(ψ) : q_last - A*log(1-ψ) + A*log(1-ψ_last) + dqdpsi_hybrid(ψ) = ψ <= ψ_last ? dqdpsi_real(ψ) : A/(1-ψ) + + function find_psi_for_q(target_q) + if target_q <= q_last + ψgrid = range(x_knots[1], ψ_last, length=200_000) + vals = q_of_psi_real.(ψgrid) + for i in 1:length(ψgrid)-1 + y1, y2 = vals[i]-target_q, vals[i+1]-target_q + if y1*y2 <= 0 + lo, hi = ψgrid[i], ψgrid[i+1] + flo = y1 + for _ in 1:60 + mid = (lo+hi)/2 + fmid = q_of_psi_real(mid) - target_q + sign(fmid) == sign(flo) ? (lo,flo)=(mid,fmid) : (hi=mid) + end + return (lo+hi)/2, true + end + end + return nothing, false + else + x = (1-ψ_last) * exp((q_last - target_q)/A) + x < 1e-15 && return nothing, false + return 1.0 - x, true + end + end + + kfile = joinpath("examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode_kinetic.h5") + psi_k, ne_k, Ti_k, Te_k, omE_k, chie_k, chip_k = h5open(kfile, "r") do kf + (read(kf["psi"]), read(kf["n_e"]), read(kf["T_i"]), read(kf["T_e"]), + read(kf["omega_E"]), read(kf["chi_e"]), read(kf["chi_phi"])) + end + + kp = KineticProfiles( + psi = psi_k, + n_e = ne_k, + T_e = Te_k, + T_i = Ti_k, + omega = omE_k, + omega_e = zeros(length(psi_k)), + omega_i = zeros(length(psi_k)), + ) + chie_spl = cubic_interp(psi_k, chie_k; extrap=Extrap(:extend)) + chip_spl = cubic_interp(psi_k, chip_k; extrap=Extrap(:extend)) + + chi1 = 2π * equil.psio + function omega_star_at(ψ) + n_e_val = kp.n_e(ψ); dn_e = kp.n_e(ψ; deriv=DerivOp(1)) + T_e_val = kp.T_e(ψ); dT_e = kp.T_e(ψ; deriv=DerivOp(1)) + T_i_val = kp.T_i(ψ); dT_i = kp.T_i(ψ; deriv=DerivOp(1)) + ω_star_e = (2π / chi1) * (T_e_val * dn_e / n_e_val + dT_e) + ω_star_i = -(2π / (z_i * chi1)) * (T_i_val * dn_e / n_e_val + dT_i) + return (ω_star_e, ω_star_i) + end + + R0 = equil.ro + bt_at(ψ) = Float64(equil.profiles.F_spline(ψ)) / (2π * R0) + + kept_r, kept_q, kept_psi = Float64[], Float64[], Float64[] + kept_w_dels, kept_w_visco = Float64[], Float64[] + psihigh_dels, psihigh_visco = nothing, nothing + + for m in q_scan + ψ_m, found = find_psi_for_q(Float64(m)) + !found && continue + + rs_val = surface_minor_radius(equil, ψ_m) + q_val = q_of_psi_hybrid(ψ_m) + q1_val = dqdpsi_hybrid(ψ_m) + da_dpsi = surface_da_dpsi(equil, ψ_m) + (!isfinite(da_dpsi) || da_dpsi == 0.0 || !isfinite(q1_val) || !isfinite(q_val)) && continue + + sval_r = try + r_based_shear(rs_val, q_val, q1_val, da_dpsi) + catch + continue + end + + prof = kp(ψ_m) + ω_e_use, ω_i_use = omega_star_at(ψ_m) + chi_perp, chi_tor = chie_spl(ψ_m), chip_spl(ψ_m) + (!isfinite(chi_perp) || !isfinite(chi_tor) || chi_perp<=0 || chi_tor<=0) && continue + + p = try + 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, + qval=q_val, sval_r=sval_r, bt=bt_at(ψ_m), rs=rs_val, R0=R0, + mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, + m=Int(round(m)), n=n_tor, ising=length(kept_q)+1) + catch + continue + end + + lw = slayer_layer_thickness(p) + isnan(lw.delta_s_m) && continue + + push!(kept_r, rs_val); push!(kept_q, m); push!(kept_psi, ψ_m) + push!(kept_w_dels, lw.delta_s_m); push!(kept_w_visco, lw.δ_visco) + + n_now = length(kept_psi) + if n_now >= 2 + right_prev_dels = kept_psi[n_now-1] + kept_w_dels[n_now-1]/2 + left_curr_dels = kept_psi[n_now] - kept_w_dels[n_now]/2 + right_prev_visco = kept_psi[n_now-1] + kept_w_visco[n_now-1]/2 + left_curr_visco = kept_psi[n_now] - kept_w_visco[n_now]/2 + + dels_stop = (left_curr_dels < right_prev_dels) || (right_prev_dels > 1.0) + visco_stop = (left_curr_visco < right_prev_visco) || (right_prev_visco > 1.0) + + if dels_stop && psihigh_dels === nothing + psihigh_dels = kept_psi[n_now-1] - kept_w_dels[n_now-1]/2 + end + if visco_stop && psihigh_visco === nothing + psihigh_visco = kept_psi[n_now-1] - kept_w_visco[n_now-1]/2 + end + dels_stop && visco_stop && break + end + end + + candidates = filter(!isnothing, [psihigh_dels, psihigh_visco]) + psihigh_final = isempty(candidates) ? nothing : minimum(candidates) + + print_layer_width_table(kept_q, kept_psi, kept_w_dels, "Resistive layer width (Riccati ODE)") + print_layer_width_table(kept_q, kept_psi, kept_w_visco, "Visco-resistive comparison scale") + + return psihigh_final, psihigh_dels, psihigh_visco +end + +function print_layer_width_table(kept_q, kept_psi, kept_w, label) + println("\n$label") + println(rpad("q", 6), rpad("psi", 14), rpad("left edge", 14), rpad("right edge", 14), "width") + for i in eachindex(kept_q) + left_edge = kept_psi[i] - kept_w[i]/2 + right_edge = kept_psi[i] + kept_w[i]/2 + println( + rpad(string(kept_q[i]), 6), + rpad(string(round(kept_psi[i], digits=8)), 14), + rpad(string(round(left_edge, digits=8)), 14), + rpad(string(round(right_edge, digits=8)), 14), + string(round(kept_w[i], digits=10)) + ) + end +end