Merge pull request #40 from ProjectTorreyPines/multicoil

Implement MultiCoil

src/VacuumFields.jl

@@ -22,7 +22,7 @@ const IMASelement = IMAS.pf_active__coil___element
 const IMASoutline = Union{IMAS.pf_active__coil___element___geometry__outline, NamedTuple}
 
 include("coils.jl")
-export AbstractCoil, PointCoil, ParallelogramCoil, QuadCoil, DistributedCoil, encircling_coils
+export AbstractCoil, PointCoil, ParallelogramCoil, QuadCoil, DistributedCoil, MultiCoil, encircling_coils
 
 
 include("integration.jl")

src/coils.jl

@@ -1,26 +1,33 @@
 """
-    AbstractCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real}
+    AbstractCoil
 
 Abstract coil type
 """
-abstract type AbstractCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} end
-(::Type{T})(R, Z; resistance=0.0, turns=1) where {T<:AbstractCoil} = T(R, Z, 0.0, resistance, turns)
-(::Type{T})(R, Z, current; resistance=0.0, turns=1) where {T<:AbstractCoil} = T(R, Z, current, resistance, turns)
+abstract type AbstractCoil end
 
 """
-    PointCoil{T1, T2, T3, T4} <:  AbstractCoil{T1, T2, T3, T4}
+    AbstractSingleCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real}
+
+Abstract coil type
+"""
+abstract type AbstractSingleCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractCoil end
+(::Type{T})(R, Z; resistance=0.0, turns=1) where {T<:AbstractSingleCoil} = T(R, Z, 0.0, resistance, turns)
+(::Type{T})(R, Z, current; resistance=0.0, turns=1) where {T<:AbstractSingleCoil} = T(R, Z, current, resistance, turns)
+
+"""
+    PointCoil{T1, T2, T3, T4} <:  AbstractSingleCoil{T1, T2, T3, T4}
 
 Point filament coil at scalar (R, Z) location
 """
-mutable struct PointCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractCoil{T1,T2,T3,T4}
+mutable struct PointCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractSingleCoil{T1,T2,T3,T4}
     R::T1
     Z::T1
     current::T2
     resistance::T3
     turns::T4
 end
 
-current(coil::AbstractCoil) = coil.current
+current(coil::AbstractSingleCoil) = coil.current
 
 # BCL 2/27/24
 # N.B.: Not sure about sign with turns and such
@@ -33,16 +40,32 @@ function current(coil::IMAScoil)
     end
 end
 
+@inline function set_current!(coil::AbstractSingleCoil, current::Real)
+    coil.current = current
+end
+
 function set_current!(coil::IMAScoil, current::Real)
     turns = sum(element.turns_with_sign for element in coil.element)
     ipt = current / turns
     IMAS.@ddtime(coil.current.data = ipt)
     return coil
 end
 
-resistance(coil::Union{AbstractCoil,IMAScoil}) = coil.resistance
+resistance(coil::AbstractSingleCoil) = coil.resistance
+
+function resistance(coil::IMAScoil{T}) where {T<:Real}
+    try
+        return coil.resistance
+    catch e
+        if e isa IMAS.IMASmissingDataException
+            return zero(T)
+        else
+            rethrow(e)
+        end
+    end
+end
 
-turns(coil::AbstractCoil) = coil.turns
+turns(coil::AbstractSingleCoil) = coil.turns
 # VacuumFields turns are sign-less
 turns(coil::IMAScoil) = abs(sum(element.turns_with_sign for element in coil.element))
 turns(element::IMASelement) = abs(element.turns_with_sign)
@@ -56,12 +79,12 @@ Return PointCoil, a point filament coil at scalar (R, Z) location
 PointCoil(R, Z, current=0.0; resistance=0.0, turns=1) = PointCoil(R, Z, current, resistance, turns)
 
 """
-    ParallelogramCoil{T1, T2, T3, T4} <:  AbstractCoil{T1, T2, T3, T4}
+    ParallelogramCoil{T1, T2, T3, T4} <:  AbstractSingleCoil{T1, T2, T3, T4}
 
 Parallelogram coil with the R, Z, ΔR, ΔZ, θ1, θ2 formalism (as used by EFIT, for example).
 Here θ1 and θ2 are the shear angles along the x- and y-axes, respectively, in degrees.
 """
-mutable struct ParallelogramCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractCoil{T1,T2,T3,T4}
+mutable struct ParallelogramCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractSingleCoil{T1,T2,T3,T4}
     R::T1
     Z::T1
     ΔR::T1
@@ -89,11 +112,11 @@ function area(ΔR, ΔZ, θ1, θ2)
 end
 
 """
-    QuadCoil{T1, T2, T3, T4} <:  AbstractCoil{T1, T2, T3, T4}
+    QuadCoil{T1, T2, T3, T4} <:  AbstractSingleCoil{T1, T2, T3, T4}
 
 Quadrilateral coil with counter-clockwise corners (starting from lower left) at R and Z
 """
-mutable struct QuadCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real,VT1<:AbstractVector{T1}} <: AbstractCoil{T1,T2,T3,T4}
+mutable struct QuadCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real,VT1<:AbstractVector{T1}} <: AbstractSingleCoil{T1,T2,T3,T4}
     R::VT1
     Z::VT1
     current::T2
@@ -129,6 +152,15 @@ function QuadCoil(pc::ParallelogramCoil)
     return QuadCoil(R, Z, pc.current; pc.resistance, pc.turns)
 end
 
+
+function QuadCoil(elm::IMASelement, current_per_turn::Real = 0.0, resistance_per_turn::Real = 0.0)
+    R, Z = IMAS.outline(elm)
+    @assert length(R) == length(Z) == 4
+    Nt = turns(elm)
+    return QuadCoil(R, Z, current_per_turn * Nt, resistance_per_turn * Nt, Nt)
+end
+
+
 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]
@@ -162,11 +194,11 @@ function resistance(element::IMASelement, resistivity::Real)
 end
 
 """
-    DistributedCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractCoil{T1,T2,T3,T4}
+    DistributedCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractSingleCoil{T1,T2,T3,T4}
 
 Coil consisting of distributed filaments
 """
-mutable struct DistributedCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractCoil{T1,T2,T3,T4}
+mutable struct DistributedCoil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractSingleCoil{T1,T2,T3,T4}
     R::Vector{T1}
     Z::Vector{T1}
     current::T2
@@ -227,6 +259,49 @@ function DistributedCoil(C::ParallelogramCoil; spacing::Real=0.01)
     return DistributedParallelogramCoil(C.R, C.Z, C.ΔR, C.ΔZ, C.θ1, C.θ2, C.current; spacing, C.turns)
 end
 
+
+mutable struct MultiCoil{SC<:AbstractSingleCoil} <: AbstractCoil
+    coils::Vector{SC}
+end
+
+function MultiCoil(icoil::IMAScoil)
+    total_turns = turns(icoil)
+    resistance_per_turn = resistance(icoil) / total_turns
+
+    current_per_turn = current(icoil) / total_turns
+    return MultiCoil([QuadCoil(elm, current_per_turn, resistance_per_turn) for elm in icoil.element])
+end
+
+
+function is_active(coil::IMAScoil)
+    funcs = (IMAS.index_2_name(coil.function)[f.index] for f in coil.function)
+    return (:shaping in funcs) && current(coil) != 0.0
+end
+
+function MultiCoils(dd::IMAS.dd{D}; active_only::Bool=true) where {D<:Real}
+    if active_only
+        coils = [MultiCoil(coil) for coil in filter(is_active, dd.pf_active.coil)]
+    else
+        coils = [MultiCoil(coil) for coil in dd.pf_active.coil]
+    end
+    return coils
+end
+
+current(mcoil::MultiCoil) = sum(current(coil) for coil in mcoil.coils)
+resistance(mcoil::MultiCoil) = sum(resistance(coil) for coil in mcoil.coils)
+turns(mcoil::MultiCoil) = sum(turns(coil) for coil in mcoil.coils)
+
+function set_current!(mcoil::MultiCoil, current::Real)
+    total_turns = turns(mcoil)
+    current_per_turn = current / total_turns
+    for coil in mcoil.coils
+        Nt = turns(coil)
+        set_current!(coil, current_per_turn * Nt)
+    end
+    return mcoil
+end
+
+
 """
     encircling_coils(EQfixed::MXHEquilibrium.AbstractEquilibrium, n_coils::Int)
 
@@ -381,4 +456,14 @@ end
             coil
         end
     end
+end
+
+@recipe function plot_coils(mcoil::MultiCoil)
+    for (k, coil) in enumerate(mcoil.coils)
+        @series begin
+            primary := (k == 1)
+            aspect_ratio := :equal
+            coil
+        end
+    end
 end

src/control.jl

@@ -308,7 +308,7 @@ function find_coil_currents!(
 
     # First reset current in coils to unity
     for coil in coils
-        coil.current = 1.0
+        set_current!(coil, 1.0)
     end
 
     N = length(flux_cps) + 2 * length(saddle_cps)
@@ -328,7 +328,7 @@ function find_coil_currents!(
 
     # update values of coils current
     for (k, coil) in enumerate(coils)
-        coil.current = Ic0[k]
+        set_current!(coil, Ic0[k])
     end
 
     cost = norm(A * Ic0 .- b) / norm(b)

src/fixed2free.jl

@@ -123,22 +123,22 @@ plasma_boundary_psi_w_fallback(shot::TEQUILA.Shot, args...) = MXHEquilibrium.Bou
 """
     fixed2free(
         EQfixed::MXHEquilibrium.AbstractEquilibrium,
-        coils::AbstractVector{<:AbstractCoil{T,C}},
+        coils::AbstractVector{<:AbstractCoil},
         R::AbstractVector{T},
         Z::AbstractVector{T};
         ψbound::Real=0.0,
-        kwargs...) where {T<:Real,C<:Real}
+        kwargs...) where {T<:Real}
 
 Convert the flux of a fixed-boundary equilibrium `EQfixed` to a free-boundary representation on an `(R,Z)` grid,
     using the flux from `coils` with currents satisfying given control points
 """
 function fixed2free(
     EQfixed::MXHEquilibrium.AbstractEquilibrium,
-    coils::AbstractVector{<:AbstractCoil{T,C}},
+    coils::AbstractVector{<:AbstractCoil},
     R::AbstractVector{T},
     Z::AbstractVector{T};
     ψbound::Real=0.0,
-    kwargs...) where {T<:Real,C<:Real}
+    kwargs...) where {T<:Real}
     image = Image(EQfixed)
     return fixed2free(EQfixed, image, coils, R, Z; ψbound, kwargs...)
 end
@@ -147,14 +147,14 @@ end
     fixed2free(
         EQfixed::MXHEquilibrium.AbstractEquilibrium,
         image::Image,
-        coils::AbstractVector{<:AbstractCoil{T,C}},
+        coils::AbstractVector{<:AbstractCoil},
         R::AbstractVector{T},
         Z::AbstractVector{T};
         flux_cps::Vector{<:FluxControlPoint}=FluxControlPoint{Float64}[],
         saddle_cps::Vector{<:SaddleControlPoint}=SaddleControlPoint{Float64}[],
         ψbound::Real=0.0,
         fixed_coils::Vector{<:AbstractCoil}=PointCoil{Float64,Float64}[],
-        λ_regularize::Real=0.0) where {T<:Real,C<:Real}
+        λ_regularize::Real=0.0) where {T<:Real}
 
 Convert the flux of a fixed-boundary equilibrium `EQfixed` to a free-boundary representation on an `(R,Z)` grid,
     subtracting out the flux from image currents `image` and adding in the flux from `coils` with currents
@@ -163,14 +163,14 @@ Convert the flux of a fixed-boundary equilibrium `EQfixed` to a free-boundary re
 function fixed2free(
     EQfixed::MXHEquilibrium.AbstractEquilibrium,
     image::Image,
-    coils::AbstractVector{<:AbstractCoil{T,C}},
+    coils::AbstractVector{<:AbstractCoil},
     R::AbstractVector{T},
     Z::AbstractVector{T};
     flux_cps::Vector{<:FluxControlPoint}=FluxControlPoint{Float64}[],
     saddle_cps::Vector{<:SaddleControlPoint}=SaddleControlPoint{Float64}[],
     ψbound::Real=0.0,
     fixed_coils::Vector{<:AbstractCoil}=PointCoil{Float64,Float64}[],
-    λ_regularize::Real=0.0) where {T<:Real,C<:Real}
+    λ_regularize::Real=0.0) where {T<:Real}
 
     Sb, ψb = plasma_boundary_psi_w_fallback(EQfixed)
     ψ_f2f = T[MXHEquilibrium.in_boundary(Sb, (r, z)) ? EQfixed(r, z) - ψb + ψbound : ψbound for z in Z, r in R]

src/flux.jl

@@ -9,12 +9,17 @@
 
 
 # coil flux
-@inline function _pfunc(Gfunc, coil::Union{AbstractCoil, IMAScoil}, R::Real, Z::Real;
+@inline function _pfunc(Gfunc, coil::Union{AbstractSingleCoil, IMAScoil}, R::Real, Z::Real;
                         COCOS::MXHEquilibrium.COCOS=MXHEquilibrium.cocos(11), Bp_fac::Float64=COCOS.sigma_Bp * (2π)^COCOS.exp_Bp,
                         coil_current::Real=current(coil))
     return μ₀ * Bp_fac * Gfunc(coil, R, Z) * coil_current
 end
 
+@inline function _pfunc(Gfunc, mcoil::MultiCoil, R::Real, Z::Real;
+    COCOS::MXHEquilibrium.COCOS=MXHEquilibrium.cocos(11), Bp_fac::Float64=COCOS.sigma_Bp * (2π)^COCOS.exp_Bp)
+    return μ₀ * Bp_fac * sum(Gfunc(coil, R, Z) * current(coil) for coil in mcoil.coils)
+end
+
 # image flux
 @inline function _pfunc(Gfunc, image::Image, R::Real, Z::Real)
     tid = Threads.threadid()

src/green.jl

@@ -19,10 +19,15 @@ end
 function _gfunc(Gfunc::Function, element::IMASelement, R::Real, Z::Real, scale_factor::Real=1.0; xorder::Int=default_order, yorder::Int=default_order)
     return _gfunc(Gfunc, IMAS.outline(element), R, Z, scale_factor; xorder, yorder)
 end
+
 function _gfunc(Gfunc::Function, ol::IMASoutline, R::Real, Z::Real, scale_factor::Real=1.0; xorder::Int=default_order, yorder::Int=default_order)
     return integrate((X, Y) -> Gfunc(X, Y, R, Z, scale_factor), ol; xorder, yorder) / area(ol)
 end
 
+function _gfunc(Gfunc::Function, mcoil::MultiCoil, R::Real, Z::Real, scale_factor::Real=1.0; kwargs...)
+    return sum(current(coil) * _gfunc(Gfunc, coil, R, Z, scale_factor; kwargs...) for coil in mcoil.coils) / current(mcoil)
+end
+
 # Generalized wrapper functions for all coil types
 function Green(coil::Union{AbstractCoil, IMAScoil, IMASelement}, R::Real, Z::Real, scale_factor::Real=1.0; kwargs...)
     return _gfunc(Green, coil, R, Z, scale_factor; kwargs...)

src/imas_coils.jl

@@ -1,7 +1,7 @@
 #= ==================================== =#
 #  IMAS.pf_active__coil to VacuumFields  #
 #= ==================================== =#
-mutable struct GS_IMAS_pf_active__coil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractCoil{T1,T2,T3,T4}
+mutable struct GS_IMAS_pf_active__coil{T1<:Real,T2<:Real,T3<:Real,T4<:Real} <: AbstractSingleCoil{T1,T2,T3,T4}
     imas::IMAS.pf_active__coil{T1}
     tech::IMAS.build__pf_active__technology{T1}
     time0::Float64
@@ -205,10 +205,11 @@ end
 
 @recipe function plot_coil(coil::GS_IMAS_pf_active__coil)
     @series begin
-        seriestype := :scatter
-        marker --> :circle
+        #seriestype := :scatter
+        #marker --> :circle
         label --> ""
-        [coil.r], [coil.z]
+        #[coil.r], [coil.z]
+        coil.imas
     end
 end
 

src/mutual.jl

@@ -37,6 +37,10 @@ function mutual(C1::Union{AbstractCoil, IMASelement}, C2::Union{ParallelogramCoi
     return fac * integrate(f, C2; xorder, yorder) / area(C2)
 end
 
+function mutual(C1::Union{AbstractCoil, IMASelement}, mcoil::MultiCoil; kwargs...)
+    return sum(mutual(C1, C2; kwargs...) for C2 in mcoil.coils)
+end
+
 """
     mutual(C1::Union{AbstractCoil, IMASelement}, C2::IMAScoil; xorder::Int=3, yorder::Int=3)