From e4c4b4f3f1195b5d802ac979aee368b2870fb8ef Mon Sep 17 00:00:00 2001 From: Nan Shi Date: Tue, 24 Mar 2026 15:45:01 -0700 Subject: [PATCH 1/2] Fix sonic/toroidal rotation conversion near axis and use thermal pressure MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Smooth dp/dψ correction to zero near ρ<0.2 with quadratic ramp to avoid singularity at magnetic axis. Use pressure_thermal instead of total pressure (which includes fast ions) for the diamagnetic correction. --- src/physics/neoclassical.jl | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/src/physics/neoclassical.jl b/src/physics/neoclassical.jl index 24308edf..31419396 100644 --- a/src/physics/neoclassical.jl +++ b/src/physics/neoclassical.jl @@ -492,8 +492,10 @@ function ωtor2sonic(cp1d::IMAS.core_profiles__profiles_1d{T}; ind::Int=0) where ni = cp1d.ion[ind].density zi = cp1d.ion[ind].z_ion rot_freq = cp1d.ion[ind].rotation_frequency_tor - dpi_dpsi = gradient(cp1d.grid.psi, cp1d.ion[ind].pressure) - return @. rot_freq + dpi_dpsi / (ni * zi * IMAS.mks.e) + rho = cp1d.grid.rho_tor_norm + dpi_dpsi = gradient(cp1d.grid.psi, cp1d.ion[ind].pressure_thermal) + correction = @. dpi_dpsi / (ni * zi * IMAS.mks.e) * min(rho / 0.2, 1.0)^2 + return @. rot_freq + correction end @compat public ωtor2sonic @@ -522,8 +524,10 @@ function sonic2ωtor(cp1d::IMAS.core_profiles__profiles_1d, ion::IMAS.core_profi ωExB = cp1d.rotation_frequency_tor_sonic ni = ion.density zi = ion.z_ion - dpi_dpsi = gradient(cp1d.grid.psi, ion.pressure) - return @. ωExB - dpi_dpsi / (ni * zi * IMAS.mks.e) + rho = cp1d.grid.rho_tor_norm + dpi_dpsi = gradient(cp1d.grid.psi, ion.pressure_thermal) + correction = @. dpi_dpsi / (ni * zi * IMAS.mks.e) * min(rho / 0.2, 1.0)^2 + return @. ωExB - correction end @compat public sonic2ωtor From c6a0cd146a9265fb665365532d5b2ee22bb39abf Mon Sep 17 00:00:00 2001 From: jmcclena Date: Tue, 24 Mar 2026 20:15:30 -0700 Subject: [PATCH 2/2] Update neoclassical.jl minor updates to on-axis dpi_dpsi correction factor --- src/physics/neoclassical.jl | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/physics/neoclassical.jl b/src/physics/neoclassical.jl index 31419396..1b5c8c37 100644 --- a/src/physics/neoclassical.jl +++ b/src/physics/neoclassical.jl @@ -481,20 +481,23 @@ Populates ExB rotation based on ion omega_tor assuming ωpol=0.0 Returns the updated `cp1d` with `rotation_frequency_tor_sonic` populated and all ion rotation frequencies updated. `ind` is the ion species index to use as reference (0 = auto-select first non-zero rotation) + +`rho_correction` is the near axis value where begin to taper dpi_dpsi term to zero, + (This is needed to handle large gradients from fast ion subtraction) """ -function ωtor2sonic(cp1d::IMAS.core_profiles__profiles_1d{T}; ind::Int=0) where {T<:Real} +function ωtor2sonic(cp1d::IMAS.core_profiles__profiles_1d{T}; ind::Int=0, rho_correction::Real=0.2) where {T<:Real} if ind == 0 ind = findfirst(x -> hasdata(x, :rotation_frequency_tor) && sum(abs.(x.rotation_frequency_tor)) != 0.0, cp1d.ion) if ind === nothing return zero(cp1d.grid.rho_tor_norm) end end - ni = cp1d.ion[ind].density + ni = cp1d.ion[ind].density_thermal zi = cp1d.ion[ind].z_ion rot_freq = cp1d.ion[ind].rotation_frequency_tor rho = cp1d.grid.rho_tor_norm dpi_dpsi = gradient(cp1d.grid.psi, cp1d.ion[ind].pressure_thermal) - correction = @. dpi_dpsi / (ni * zi * IMAS.mks.e) * min(rho / 0.2, 1.0)^2 + correction = @. dpi_dpsi / (ni * zi * IMAS.mks.e) * min(rho / rho_correction, 1.0)^2 return @. rot_freq + correction end @@ -518,15 +521,18 @@ push!(document[Symbol("Physics neoclassical")], :ωtor2sonic!) Returns ion ω_tor based on ExB rotation assuming ω_pol=0.0 +`rho_correction` is the near axis value where begin to taper dpi_dpsi term to zero, + (This is needed to handle large gradients from fast ion subtraction) + Uses the existing `rotation_frequency_tor_sonic` field to compute individual ion toroidal rotation frequencies """ -function sonic2ωtor(cp1d::IMAS.core_profiles__profiles_1d, ion::IMAS.core_profiles__profiles_1d___ion) +function sonic2ωtor(cp1d::IMAS.core_profiles__profiles_1d, ion::IMAS.core_profiles__profiles_1d___ion; rho_correction::Real=0.2) ωExB = cp1d.rotation_frequency_tor_sonic - ni = ion.density + ni = ion.density_thermal zi = ion.z_ion rho = cp1d.grid.rho_tor_norm dpi_dpsi = gradient(cp1d.grid.psi, ion.pressure_thermal) - correction = @. dpi_dpsi / (ni * zi * IMAS.mks.e) * min(rho / 0.2, 1.0)^2 + correction = @. dpi_dpsi / (ni * zi * IMAS.mks.e) * min(rho / rho_correction, 1.0)^2 return @. ωExB - correction end