diff --git a/src/physics/pedestal.jl b/src/physics/pedestal.jl index 395da686..78508a69 100644 --- a/src/physics/pedestal.jl +++ b/src/physics/pedestal.jl @@ -52,8 +52,10 @@ function blend_core_edge_Hmode( p_targets = interp1d(rho, profile).(rho_targets) # figure out expin and expout such that the Z's of Hmode_profiles match the z_targets from transport - x_guess = [1.0, 1.0] - res = Optim.optimize(x -> cost_find_EPED_exps(x, ped_height, ped_width, rho, profile, p_targets, z_targets, rho_targets), x_guess, Optim.NelderMead()) + expin = 1.0 + expout = 1.0 + x_guess = [expin, expout] + res = Optim.optimize(x -> cost_find_EPED_exps(x, ped_height, ped_width, rho, profile, p_targets, z_targets, rho_targets), x_guess, Optim.NelderMead(), Optim.Options(;g_tol=1E-6)) expin = abs(res.minimizer[1]) expout = abs(res.minimizer[2]) @@ -75,7 +77,7 @@ push!(document[Symbol("Physics pedestal")], :blend_core_edge_Hmode) expout::Real ) -Blends the core and pedestal for given profile to match ped_height, ped_width using nml_bound as blending boundary +Blends the core and pedestal for given profile to match ped_height, ped_width using nml_bound and ped_bound as blending boundaries """ function blend_core_edge_EPED( profile::AbstractVector{<:Real}, @@ -86,20 +88,43 @@ function blend_core_edge_EPED( ped_bound::Real, expin::Real, expout::Real +) + # H-mode profile used for pedestal + # NOTE: Note that we do not provide a core value as this causes the pedestal solution to depend on the core solution + # which breaks finding self-consistent core-pedestal solution through an optimizer + profile_ped = ped_height_at_09(rho, Hmode_profiles(profile[end], ped_height, length(rho), expin, expout, ped_width), ped_height) + return blend_core_edge(profile, profile_ped, rho, nml_bound, ped_bound) +end + +@compat public blend_core_edge_EPED +push!(document[Symbol("Physics pedestal")], :blend_core_edge_EPED) + +""" + blend_core_edge( + profile::AbstractVector{<:Real}, + profile_ped::AbstractVector{<:Real}, + rho::AbstractVector{<:Real}, + nml_bound::Real, + ped_bound::Real) + +Blends two profiles core and edge using nml_bound and ped_bound as blending boundaries +""" +function blend_core_edge( + profile::AbstractVector{<:Real}, + profile_ped::AbstractVector{<:Real}, + rho::AbstractVector{<:Real}, + nml_bound::Real, + ped_bound::Real ) @assert rho[end] == 1.0 @assert nml_bound <= ped_bound "Unable to blend the core-pedestal because the nml_bound $nml_bound > ped_bound top $ped_bound" + @assert length(profile) == length(profile_ped) == length(rho) iped = argmin(abs.(rho .- ped_bound)) inml = argmin(abs.(rho .- nml_bound)) z_profile = -calc_z(rho, profile, :backward) z_nml = z_profile[inml] - # H-mode profile used for pedestal - # NOTE: Note that we do not provide a core value as this causes the pedestal solution to depend on the core solution - # which breaks finding self-consistent core-pedestal solution through an optimizer - profile_ped = Hmode_profiles(profile[end], ped_height, length(rho), expin, expout, ped_width) - # linear z between nml and pedestal if nml_bound < ped_bound z_profile_ped = -calc_z(rho, profile_ped, :backward) @@ -118,9 +143,6 @@ function blend_core_edge_EPED( return profile_new end -@compat public blend_core_edge_EPED -push!(document[Symbol("Physics pedestal")], :blend_core_edge_EPED) - """ blend_core_edge_Lmode( profile::AbstractVector{<:Real}, @@ -250,3 +272,15 @@ end @compat public pedestal_finder push!(document[Symbol("Physics pedestal")], :pedestal_finder) + +""" + ped_height_at_09(rho::AbstractVector{T}, profile::AbstractVector{T}, height09::T) where {T<:Real} + +Scale density profile so that value at rho=0.9 is height09 +""" +function ped_height_at_09(rho::AbstractVector{T}, profile::AbstractVector{T}, ped_height::T) where {T<:Real} + return profile / interp1d(rho, profile).(0.9) * ped_height +end + +@compat public ped_height_at_09 +push!(document[Symbol("Physics pedestal")], :ped_height_at_09)