merging dev into master

ProjectTorreyPines · Nov 10, 2024 · bb09f08 · bb09f08
2 parents 931e614 + ca8181c
commit bb09f08
Show file tree

Hide file tree

Showing 4 changed files with 52 additions and 48 deletions.
diff --git 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)
@@ -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
 
@@ -558,54 +560,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

diff --git 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
 
 

diff --git 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
@@ -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)