Skip to content
Open
Changes from all commits
Commits
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
18 changes: 9 additions & 9 deletions src/physics/transport.jl
Original file line number Diff line number Diff line change
Expand Up @@ -61,11 +61,11 @@ If `rho_ped > transport_grid[end]` then rotation shear is linearly interpolated

If `rho_ped < transport_grid[end]` then boundary condition is at `transport_grid[end]`
"""
function profile_from_rotation_shear_transport(
function profile_from_grad(
profile_old::AbstractVector{<:Real},
rho::AbstractVector{<:Real},
transport_grid::AbstractVector{<:Real},
rotation_shear_transport_grid::AbstractVector{<:Real},
grad_transport_grid::AbstractVector{<:Real},
rho_ped::Real=0.0)

transport_indices = [IMAS.argmin_abs(rho, rho_x) for rho_x in transport_grid]
Expand All @@ -74,17 +74,17 @@ function profile_from_rotation_shear_transport(

if index_ped > index_last
# If rho_ped is beyond transport_grid[end], extend interpolation
dw_dr_old = gradient(rho, profile_old; method=:backward)
grad_old = grad(rho, profile_old; method=:backward)
transport_indices = vcat(1, transport_indices, index_ped)
rotation_shear_transport_grid = vcat(0.0, rotation_shear_transport_grid, dw_dr_old[index_ped])
grad_transport_grid = vcat(0.0, grad_transport_grid, grad_old[index_ped])
else
# If rho_ped is within transport_grid, use transport_grid[end] as boundary
transport_indices = vcat(1, transport_indices)
rotation_shear_transport_grid = vcat(0.0, rotation_shear_transport_grid)
grad_transport_grid = vcat(0.0, grad_transport_grid)
end

# Interpolate rotation shear to the integration range
dw_dr = IMAS.interp1d(transport_indices, rotation_shear_transport_grid).(1:index_last)
grad = IMAS.interp1d(transport_indices, grad_transport_grid).(1:index_last)

# Set up the profile
profile_new = similar(profile_old)
Expand All @@ -95,9 +95,9 @@ function profile_from_rotation_shear_transport(
profile_new[index_last] = profile_old[index_last]
for i in (index_last-1):-1:1
# Trapezoidal integration: negative because we integrate inward
dw_avg = (dw_dr[i] + dw_dr[i+1]) / 2.0
grad_avg = (grad[i] + grad[i+1]) / 2.0
dr = rho[i+1] - rho[i]
profile_new[i] = profile_new[i+1] - dw_avg * dr
profile_new[i] = profile_new[i+1] - grad_avg * dr
end

return profile_new
Expand Down Expand Up @@ -281,4 +281,4 @@ end


@compat public total_fluxes
push!(document[Symbol("Physics transport")], :total_fluxes)
push!(document[Symbol("Physics transport")], :total_fluxes)