Skip to content

Fix sonic/toroidal rotation conversion near axis and use thermal pressure #296

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jmcclena merged 2 commits into from
Mar 25, 2026
Merged

Fix sonic/toroidal rotation conversion near axis and use thermal pressure #296

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
e4c4b4f
Fix sonic/toroidal rotation conversion near axis and use thermal pres…
nanshi1177 Mar 24, 2026
c6a0cd1
Update neoclassical.jl
jmcclena Mar 25, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
26 changes: 18 additions & 8 deletions src/physics/neoclassical.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -481,19 +481,24 @@ 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
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 / rho_correction, 1.0)^2
return @. rot_freq + correction
end

@compat public ωtor2sonic
Expand All @@ -516,14 +521,19 @@ 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
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 / rho_correction, 1.0)^2
return @. ωExB - correction
end

@compat public sonic2ωtor
Expand Down
Loading