From 59a0eb5031fb21c19590044b141aa37dd2c9c765 Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Fri, 8 Nov 2024 14:02:36 -0800 Subject: [PATCH 1/6] fix issue with negative resistance --- src/coils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coils.jl b/src/coils.jl index e755e66..2a112f5 100644 --- a/src/coils.jl +++ b/src/coils.jl @@ -232,7 +232,7 @@ area(ol::IMASoutline) = area(ol.r, ol.z) # compute the resistance given a resistitivity function resistance(coil::Union{ParallelogramCoil,QuadCoil,IMASelement}, resistivity::Real) - return 2π * turns(coil)^2 * resistivity / integrate((R, Z) -> 1.0 / R, coil) + return 2π * turns(coil)^2 * resistivity / abs(integrate((R, Z) -> 1.0 / R, coil)) end function resistance(coil::IMAScoil, resistivity::Real, element_connection::Symbol=:series) From 9affe59e152d90350ae860449935c186cabce8fc Mon Sep 17 00:00:00 2001 From: Orso Meneghini Date: Fri, 8 Nov 2024 14:03:18 -0800 Subject: [PATCH 2/6] Handle plots in cases with zero currents --- src/coils.jl | 86 +++++++++++++++++++++++++++------------------------- 1 file changed, 44 insertions(+), 42 deletions(-) diff --git a/src/coils.jl b/src/coils.jl index 2a112f5..b5f099b 100644 --- a/src/coils.jl +++ b/src/coils.jl @@ -558,54 +558,56 @@ end end max_current(coil::AbstractCoil) = abs(current(coil)) -max_current(mcoil::MultiCoil) = maximum(abs, current(coil) for coil in mcoil.coils) -max_current(icoil::IMAScoil) = abs(IMAS.@ddtime(icoil.current.data[index]) * maximum(turns(elm) for elm in icoil.element)) - -@recipe function plot_coils(coils::AbstractVector{<:Union{AbstractCoil,IMAScoil}}; color_by=:current, cname=:diverging) - @assert color_by in (nothing, :current) - if color_by === :current - alpha --> nothing - currents = [current(coil) for coil in coils] - CURRENT = maximum(max_current(coil) for coil in coils) - if CURRENT > 1e6 - currents = currents .* 1e-6 - CURRENT = CURRENT .* 1e-6 - c_unit = "MA" - elseif CURRENT > 1e3 - currents = currents .* 1e-3 - CURRENT = CURRENT .* 1e-3 - c_unit = "kA" - else - c_unit = "A" - end +max_current(mcoil::MultiCoil) = isempty(mcoil.coils) ? 0.0 : maximum(abs, current(coil) for coil in mcoil.coils) +max_current(icoil::IMAScoil) = ismissing(icoil.current, :data) ? 0.0 : (abs(IMAS.@ddtime(icoil.current.data) * maximum(turns(elm) for elm in icoil.element))) + +@recipe function plot_coils(coils::AbstractVector{<:AbstractCoil}; color_by=:current, cname=:diverging) + if !isempty(coils) + @assert color_by in (nothing, :current) + if color_by === :current + alpha --> nothing + currents = [current(coil) for coil in coils] + CURRENT = maximum(max_current(coil) for coil in coils) + if CURRENT > 1e6 + currents = currents .* 1e-6 + CURRENT = CURRENT .* 1e-6 + c_unit = "MA" + elseif CURRENT > 1e3 + currents = currents .* 1e-3 + CURRENT = CURRENT .* 1e-3 + c_unit = "kA" + else + c_unit = "A" + end - @series begin - colorbar_title := " \nPF currents [$c_unit]" - seriestype --> :scatter - color --> cname - clim --> (-CURRENT, CURRENT) - marker_z --> [-CURRENT, CURRENT] - [(NaN, NaN), (NaN, NaN)] + @series begin + colorbar_title := " \nPF currents [$c_unit]" + seriestype --> :scatter + color --> cname + clim --> (-CURRENT, CURRENT) + marker_z --> [-CURRENT, CURRENT] + [(NaN, NaN), (NaN, NaN)] + end end - end - for (k, coil) in enumerate(coils) - @series begin - if color_by === :current - colorbar_title := " \nPF currents [$c_unit]" - colorbar_entry := false - if CURRENT == 0.0 - color --> :black + for (k, coil) in enumerate(coils) + @series begin + if color_by === :current + colorbar_title := " \nPF currents [$c_unit]" + colorbar_entry := false + if CURRENT == 0.0 + color --> :black + else + current_color_index = (currents[k] + CURRENT) / (2 * CURRENT) + color --> PlotUtils.cgrad(cname)[current_color_index] + end else - current_color_index = (currents[k] + CURRENT) / (2 * CURRENT) - color --> PlotUtils.cgrad(cname)[current_color_index] + primary := (k == 1) end - else - primary := (k == 1) + aspect_ratio := :equal + label --> "" + coil end - aspect_ratio := :equal - label --> "" - coil end end end From d8358105579165e10162026be401ee1b1ea475a0 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Fri, 8 Nov 2024 16:31:31 -0800 Subject: [PATCH 3/6] Revert "fix issue with negative resistance" This reverts commit 59a0eb5031fb21c19590044b141aa37dd2c9c765. --- src/coils.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/coils.jl b/src/coils.jl index b5f099b..20e801a 100644 --- a/src/coils.jl +++ b/src/coils.jl @@ -232,7 +232,7 @@ area(ol::IMASoutline) = area(ol.r, ol.z) # compute the resistance given a resistitivity function resistance(coil::Union{ParallelogramCoil,QuadCoil,IMASelement}, resistivity::Real) - return 2π * turns(coil)^2 * resistivity / abs(integrate((R, Z) -> 1.0 / R, coil)) + return 2π * turns(coil)^2 * resistivity / integrate((R, Z) -> 1.0 / R, coil) end function resistance(coil::IMAScoil, resistivity::Real, element_connection::Symbol=:series) From 7ba9e6fce1a97fe28bd8e37f510e89797782ebd9 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Fri, 8 Nov 2024 16:33:09 -0800 Subject: [PATCH 4/6] Fix signs of Jacobian and area Works for clockwise or ccw points now --- src/coils.jl | 2 +- src/integration.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/coils.jl b/src/coils.jl index 20e801a..68a0b48 100644 --- a/src/coils.jl +++ b/src/coils.jl @@ -222,7 +222,7 @@ function area(R::AbstractVector{<:Real}, Z::AbstractVector{<:Real}) @assert length(R) == length(Z) == 4 A = R[1] * Z[2] + R[2] * Z[3] + R[3] * Z[4] + R[4] * Z[1] A -= R[2] * Z[1] + R[3] * Z[2] + R[4] * Z[3] + R[1] * Z[4] - return 0.5 * A + return 0.5 * abs(A) end area(C::QuadCoil) = area(C.R, C.Z) diff --git a/src/integration.jl b/src/integration.jl index b51024f..ad996df 100644 --- a/src/integration.jl +++ b/src/integration.jl @@ -140,7 +140,7 @@ function Jacobian(x::Real, y::Real, R::AbstractVector{<:Real}, Z::AbstractVector dZdy = a + b * x # All these need to be divided by 4 - return 0.0625 * (dRdx * dZdy - dRdy * dZdx) + return 0.0625 * abs(dRdx * dZdy - dRdy * dZdx) end From cf6519a059ec54f13952012da18a18de20261206 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Fri, 8 Nov 2024 16:34:30 -0800 Subject: [PATCH 5/6] Use current per turn in vertical stability calculations --- src/mutual.jl | 6 +++--- test/runtests.jl | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/src/mutual.jl b/src/mutual.jl index 62cfca2..9d099b1 100644 --- a/src/mutual.jl +++ b/src/mutual.jl @@ -219,7 +219,7 @@ First introduced in A. Portone, Nucl. Fusion 45 (2005) 926–932. https://doi.or """ function stability_margin(image::Image, coils::Vector{<:Union{AbstractCoil, IMAScoil}}, Ip::Real; order::Int=default_order) b = [Ip * dM_dZ(image, C, Ip) for C in coils] - K = Ip * sum(current(C) * d2M_dZ2(image, C, Ip) for C in coils) + K = Ip * sum(current(C) / turns(C) * d2M_dZ2(image, C, Ip) for C in coils) M = zeros(length(coils), length(coils)) for j in eachindex(coils) for k in eachindex(coils) @@ -257,7 +257,7 @@ This is the massless approximation and only use the passive conductors for compu """ function normalized_growth_rate(image::Image, coils::Vector{<:Union{AbstractCoil, IMAScoil}}, Ip::Real; order::Int=default_order) b = [Ip * dM_dZ(image, C, Ip) for C in coils] - K = Ip * sum(current(C) * d2M_dZ2(image, C, Ip) for C in coils) + K = Ip * sum(current(C) / turns(C) * d2M_dZ2(image, C, Ip) for C in coils) M = zeros(length(coils), length(coils)) for j in eachindex(coils) for k in eachindex(coils) @@ -277,7 +277,7 @@ function normalized_growth_rate(image::Image, coils::Vector{<:Union{AbstractCoil # reuse b vector for resistances for j in eachindex(coils) - b[j] = coils[j].resistance + b[j] = resistance(coils[j]) end R = Diagonal(b) diff --git a/test/runtests.jl b/test/runtests.jl index 9ee8ef1..983611c 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -74,7 +74,7 @@ end EQ = efit(transform_cocos(geqdsk, cc0, 11), 11) vars = matread((@__DIR__) * "/equilibria/d3d_vs_eq3_.mat") - matcoils = [ParallelogramCoil(load_par(vars["fcdata"][:,k])..., vars["Ic"][2+k]; resistance=vars["Rc"][k], turns=Int(vars["fcnturn"][k])) for k in eachindex(vars["fcnturn"])] + matcoils = [ParallelogramCoil(load_par(vars["fcdata"][:,k])..., vars["Ic"][2+k] * Int(vars["fcnturn"][k]); resistance=vars["Rc"][k], turns=Int(vars["fcnturn"][k])) for k in eachindex(vars["fcnturn"])] append!(matcoils, ParallelogramCoil(load_par(vars["vvdata"][:,k])...; resistance=vars["Rv"][k]) for k in eachindex(vars["Rv"])) @test isapprox(stability_margin(EQ, matcoils), vars["stability"]["m_s_noe"]; rtol=5e-2) From ca8181c950cfcd18f46ea7aa0cbbecd1d1ce0cb7 Mon Sep 17 00:00:00 2001 From: Brendan Lyons Date: Fri, 8 Nov 2024 16:35:05 -0800 Subject: [PATCH 6/6] Fix MultiCoils if pf_active or pf_passive is empty --- src/coils.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/coils.jl b/src/coils.jl index 68a0b48..b986bda 100644 --- a/src/coils.jl +++ b/src/coils.jl @@ -385,6 +385,8 @@ function MultiCoils(dd::IMAS.dd{D}; load_pf_active::Bool=true, active_only::Bool passive_coils = MultiCoils(dd.pf_passive) !load_pf_active && return passive_coils end + isempty(active_coils) && return passive_coils + isempty(passive_coils) && return active_coils return vcat(active_coils, passive_coils) end