added set_resistive_width_based_psihigh helper for resistive-layer-width-based psihigh truncation

Author: amlakbekaluachule

src/Tearing/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl

@@ -0,0 +1,137 @@
+"""
+    set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0)
+
+Given a completed equilibrium/stability run `r`, computes the resistive-layer
+based truncation boundary: finds rational surfaces via the iota spline,
+computes each surface's real resistive layer width (SLAYER Riccati del_s
+solve) and its visco-resistive comparison scale, and returns the more
+conservative (smaller) psihigh at which neighboring surfaces first overlap
+or a surface's own layer spans past psi=1.
+
+Returns `(psihigh, psihigh_dels, psihigh_visco)`.
+"""
+function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0)
+    r_of_psi = r.equil.geometry.avg_r_spline
+    drdpsi = ψ -> r_of_psi(ψ; deriv=DerivOp(1))
+
+    x_knots = r.equil.profiles.xs
+    q_knots = r.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_path = r.ctrl.kinetic_file
+    kf = h5open(kfile_path, "r")
+    psi_k = read(kf["psi"]); ne_k = read(kf["n_e"]); Ti_k = read(kf["T_i"])
+    Te_k = read(kf["T_e"]); omE_k = read(kf["omega_E"])
+    chie_k = read(kf["chi_e"]); chip_k = read(kf["chi_phi"])
+    close(kf)
+
+    mk_spline(y) = cubic_interp(psi_k, y; extrap=Extrap(:extend))
+    ne_spl, Ti_spl, Te_spl = mk_spline(ne_k), mk_spline(Ti_k), mk_spline(Te_k)
+    omE_spl, chie_spl, chip_spl = mk_spline(omE_k), mk_spline(chie_k), mk_spline(chip_k)
+    omega_star(ne_val, dne_dr, T_eV, B, Z) = -(T_eV*E_CHG)/(Z*E_CHG*B) * (dne_dr/ne_val)
+
+    bt = r.equil.params.bt0
+    R0 = r.equil.params.rmean
+
+    kept_psi = Float64[]; kept_w_dels = Float64[]; kept_w_visco = Float64[]
+    psihigh_dels, psihigh_visco = nothing, nothing
+
+    for m in q_scan
+        ψ_m, found = find_psi_for_q(Float64(m))
+        !found && continue
+
+        rs_val = r_of_psi(ψ_m); q_val = q_of_psi_hybrid(ψ_m); q1_val = dqdpsi_hybrid(ψ_m)
+        dr_val_geo = drdpsi(ψ_m)
+        (!isfinite(dr_val_geo) || dr_val_geo == 0.0 || !isfinite(q1_val) || !isfinite(q_val)) && continue
+
+        sval_r = try
+            r_based_shear(rs_val, q_val, q1_val, dr_val_geo)
+        catch
+            continue
+        end
+
+        ne_val = ne_spl(ψ_m); Ti_val = Ti_spl(ψ_m); Te_val = Te_spl(ψ_m)
+        (!isfinite(ne_val) || !isfinite(Ti_val) || !isfinite(Te_val) || ne_val<=0 || Ti_val<=0 || Te_val<=0) && continue
+
+        dne_dr = ne_spl(ψ_m; deriv=DerivOp(1)) / dr_val_geo
+        omega_e = omega_star(ne_val, dne_dr, Te_val, bt, -1.0)
+        omega_i = omega_star(ne_val, dne_dr, Ti_val, bt, 1.0)
+        omega = omE_spl(ψ_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=ne_val, t_e=Te_val, t_i=Ti_val,
+                omega=omega, omega_e=omega_e, omega_i=omega_i,
+                qval=q_val, sval_r=sval_r, bt=bt, 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_psi)+1)
+        catch
+            continue
+        end
+
+        lw = slayer_layer_thickness(p)
+        isnan(lw.delta_s_m) && continue
+
+        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)
+
+    return psihigh_final, psihigh_dels, psihigh_visco
+end