Merge pull request #45 from ProjectTorreyPines/series_circuit

Create AbstractCircuit and SeriesCircuit types

src/VacuumFields.jl

@@ -25,6 +25,8 @@ const IMASoutline = Union{IMAS.pf_active__coil___element___geometry__outline, Na
 include("coils.jl")
 export AbstractCoil, PointCoil, ParallelogramCoil, QuadCoil, DistributedCoil, MultiCoil, encircling_coils
 
+include("circuit.jl")
+export AbstractCircuit, SeriesCircuit, update_coil_currents!
 
 include("integration.jl")
 

src/circuit.jl

@@ -0,0 +1,31 @@
+const CoilVectorType = AbstractVector{<:Union{VacuumFields.AbstractCoil, IMAS.pf_active__coil, IMAS.pf_active__coil___element}}
+
+abstract type AbstractCircuit end
+
+mutable struct SeriesCircuit{CV<:CoilVectorType, T<:Real} <: AbstractCircuit
+    coils::CV
+    current_per_turn::T
+    signs::Vector{Int}
+    function SeriesCircuit(coils::CV, current_per_turn::T, signs; copy_coils::Bool=true) where {CV<:CoilVectorType, T<:Real}
+        @assert length(coils) == length(signs)
+        new_coils = copy_coils ? deepcopy(coils) : coils
+        update_coil_currents!(new_coils, current_per_turn, signs)
+        return new{CV, T}(new_coils, current_per_turn, signs)
+    end
+end
+
+function update_coil_currents!(series::SeriesCircuit, current_per_turn::Real=series.current_per_turn)
+    series.current_per_turn = current_per_turn
+    update_coil_currents!(series.coils, current_per_turn, series.signs)
+    return series
+end
+
+function update_coil_currents!(coils::CoilVectorType, current_per_turn::Real, signs::Vector{Int})
+    for (k, coil) in enumerate(coils)
+        current = signs[k] * turns(coil) * current_per_turn
+        set_current!(coil, current)
+    end
+    return coils
+end
+
+resistance(series::SeriesCircuit) = sum(resistance(coil) for coil in series.coils)

src/flux.jl

@@ -25,4 +25,10 @@ end
     tid = Threads.threadid()
     image._Vb[tid] .= image.dψdn_R .* Gfunc.(image.Rb, image.Zb, R, Z)
     return -trapz(image.Lb, image._Vb[tid])
-end
+end
+
+# series flux
+@inline ψ(series::SeriesCircuit, R::Real, Z::Real; kwargs...)       = sum(ψ(coil, R, Z; kwargs...) for coil in series.coils)
+@inline dψ_dR(series::SeriesCircuit, R::Real, Z::Real; kwargs...)   = sum(dψ_dR(coil, R, Z; kwargs...) for coil in series.coils)
+@inline dψ_dZ(series::SeriesCircuit, R::Real, Z::Real; kwargs...)   = sum(dψ_dZ(coil, R, Z; kwargs...) for coil in series.coils)
+@inline d2ψ_dZ2(series::SeriesCircuit, R::Real, Z::Real; kwargs...) = sum(d2ψ_dZ2(coil, R, Z; kwargs...) for coil in series.coils)

src/green.jl

@@ -49,6 +49,33 @@ function d2G_dZ2(coil::Union{AbstractCoil, IMAScoil, IMASelement}, R::Real, Z::R
     return _gfunc(d2G_dZ2, coil, R, Z, scale_factor; kwargs...)
 end
 
+# Function for circuits
+@inline function _gfunc(Pfunc, series::SeriesCircuit, R::Real, Z::Real, scale_factor::Real=1.0;
+    COCOS::MXHEquilibrium.COCOS=MXHEquilibrium.cocos(11), Bp_fac::Float64=COCOS.sigma_Bp * (2π)^COCOS.exp_Bp)
+    Ic = series.current_per_turn
+    try
+        update_coil_currents!(series, 1.0)
+        return scale_factor * Pfunc(series, R, Z; COCOS, Bp_fac) / (μ₀ * Bp_fac)
+    finally
+        update_coil_currents!(series, Ic)
+    end
+end
+
+function Green(series::SeriesCircuit, R::Real, Z::Real, scale_factor::Real=1.0; kwargs...)
+    return _gfunc(ψ, series, R, Z, scale_factor; kwargs...)
+end
+
+function dG_dR(series::SeriesCircuit, R::Real, Z::Real, scale_factor::Real=1.0; kwargs...)
+    return _gfunc(dψ_dR, series, R, Z, scale_factor; kwargs...)
+end
+
+function dG_dZ(series::SeriesCircuit, R::Real, Z::Real, scale_factor::Real=1.0; kwargs...)
+    return _gfunc(dψ_dZ, series, R, Z, scale_factor; kwargs...)
+end
+
+function d2G_dZ2(series::SeriesCircuit, R::Real, Z::Real, scale_factor::Real=1.0; kwargs...)
+    return _gfunc(d2ψ_dZ2, series, R, Z, scale_factor; kwargs...)
+end
 
 # Point-to-point Green's functions
 

src/mutual.jl

@@ -64,6 +64,18 @@ function mutual(C1::IMAScoil, C2::Union{AbstractCoil, IMAScoil, IMASelement}; xo
 end
 
 
+# Circuit
+
+function mutual(C1::Union{AbstractCoil, IMASelement}, SC2::SeriesCircuit; kwargs...)
+    return sum(SC2.signs[k] * mutual(C1, C2; kwargs...) for (k, C2) in enumerate(SC2.coils))
+end
+
+function mutual(SC1::SeriesCircuit, C2::Union{AbstractCoil, IMASelement, SeriesCircuit}; kwargs...)
+    return sum(SC1.signs[k] * mutual(C1, C2; kwargs...) for (k, C1) in enumerate(SC1.coils))
+end
+
+
+
 # Plasma-coil mutuals
 
 # Shifting plasma up by δZ is the same as shifting the coil down by δZ