Merge pull request #15 from ProjectTorreyPines/cold_start

Implement cold start from IMAS dd and store fixed_coils in Canvas and uses those by default

Project.toml

@@ -22,7 +22,7 @@ VacuumFields = "9d9223b5-c5da-4bf4-abee-2a8bb6775a49"
 DataInterpolations = "6.0.0"
 FFMPEG = "0.4"
 HypergeometricFunctions = "0.3"
-IMAS = "2.2"
+IMAS = "2.2.1"
 IMASutils = "1.2.1"
 Interpolations = "0.15"
 LoopVectorization = "0.12.171"

src/canvas.jl

@@ -30,6 +30,7 @@ mutable struct Canvas{T<:Real, VC<:CoilVectorType, I<:AbstractInterpolation, C1<
     _rs_circuit::C2
     _Ψ_at_coils::Vector{T}
     _tmp_Ncoils::Vector{T}
+    _fixed_coils::Vector{Int}
     _mutuals::Matrix{T}
     _mutuals_LU::LU{T, Matrix{T}, Vector{Int}}
     _a::Vector{T}
@@ -51,11 +52,7 @@ mutable struct Canvas{T<:Real, VC<:CoilVectorType, I<:AbstractInterpolation, C1<
     _z_cache::Vector{T}
 end
 
-function Canvas(dd::IMAS.dd{T}, Nr::Int, Nz::Int=Nr; load_pf_active=true, load_pf_passive=true) where {T<:Real}
-
-    wall_r, wall_z = IMAS.first_wall(dd.wall)
-    wall_r, wall_z = collect(wall_r), collect(wall_z)
-    #Rs, Zs = range(extrema(wall_r)..., Nr), range(extrema(wall_z)..., Nz)
+function Canvas(dd::IMAS.dd{T}, Nr::Int, Nz::Int=Nr; load_pf_active::Bool=true, load_pf_passive::Bool=true) where {T<:Real}
 
     eqt = dd.equilibrium.time_slice[]
     eqt2d = IMAS.findfirst(:rectangular, eqt.profiles_2d)
@@ -72,24 +69,51 @@ function Canvas(dd::IMAS.dd{T}, Nr::Int, Nz::Int=Nr; load_pf_active=true, load_p
         Ψ = [PSI_interpolant(r, z) for r in Rs, z in Zs]
     end
 
+    canvas = Canvas(dd, Rs, Zs, Ψ; load_pf_active, load_pf_passive)
+
+    update_bounds!(canvas)
+    trace_surfaces!(canvas)
+    gridded_Jtor!(canvas)
+    return canvas
+end
+
+function Canvas(dd::IMAS.dd{T}, Rs::StepRangeLen, Zs::StepRangeLen,
+                Ψ::Matrix{T}=zeros(T, length(Rs), length(Zs));
+                load_pf_active=true, load_pf_passive=true) where {T<:Real}
+
+    wall_r, wall_z = IMAS.first_wall(dd.wall)
+    wall_r, wall_z = collect(wall_r), collect(wall_z)
+
+    eqt = dd.equilibrium.time_slice[]
     boundary = IMAS.closed_polygon(eqt.boundary.outline.r, eqt.boundary.outline.z)
 
     # define current
     Ip = eqt.global_quantities.ip
 
     # define coils
     coils = VacuumFields.MultiCoils(dd; load_pf_active, load_pf_passive)
+    fixed_coils = Int[]
+    if load_pf_active
+        kpassive0 = length(dd.pf_active.coil)
+        for (k, coil) in enumerate(dd.pf_active.coil)
+            if :shaping ∉ (IMAS.index_2_name(coil.function)[f.index] for f in coil.function)
+                push!(fixed_coils, k)
+            end
+        end
+    else
+        kpassive0 = 0
+    end
+    if load_pf_passive
+        fixed_coils = vcat(fixed_coils, kpassive0 .+ eachindex(dd.pf_passive.loop))
+    end
 
-    Ns = length(eqt.profiles_1d.psi)
-    surfaces = Vector{IMAS.SimpleSurface{T}}(undef, Ns)
+    Nsurfaces = !ismissing(eqt.profiles_1d, :psi) ? length(eqt.profiles_1d.psi) : 129
+    surfaces = Vector{IMAS.SimpleSurface{T}}(undef, Nsurfaces)
 
-    canvas = Canvas(Rs, Zs, Ψ, Ip, coils, wall_r, wall_z, collect(boundary.r), collect(boundary.z), surfaces)
+    canvas = Canvas(Rs, Zs, Ψ, Ip, coils, wall_r, wall_z, collect(boundary.r), collect(boundary.z), surfaces; fixed_coils)
 
     set_Ψvac!(canvas)
     canvas._Ψpl .= canvas.Ψ - canvas._Ψvac
-    update_bounds!(canvas)
-    trace_surfaces!(canvas)
-    gridded_Jtor!(canvas)
 
     return canvas
 end
@@ -101,11 +125,12 @@ function Canvas(Rs::AbstractRange{T},
                 Rw::Vector{T},
                 Zw::Vector{T},
                 Rb_target::Vector{T},
-                Zb_target::Vector{T}) where {T<:Real}
+                Zb_target::Vector{T};
+                fixed_coils::Vector{Int}=Int[]) where {T<:Real}
     Nr, Nz = length(Rs), length(Zs)
     Ψ = zeros(T, Nr, Nz)
     surfaces = Vector{IMAS.SimpleSurface{T}}(undef, Nr - 1)
-    return Canvas(Rs, Zs, Ψ, Ip, coils, Rw, Zw, Rb_target, Zb_target, surfaces)
+    return Canvas(Rs, Zs, Ψ, Ip, coils, Rw, Zw, Rb_target, Zb_target, surfaces; fixed_coils)
 end
 
 function Canvas(Rs::AbstractRange{T},
@@ -117,7 +142,8 @@ function Canvas(Rs::AbstractRange{T},
                 Zw::Vector{T},
                 Rb_target::Vector{T},
                 Zb_target::Vector{T},
-                surfaces::Vector{<:IMAS.SimpleSurface}) where {T<:Real}
+                surfaces::Vector{<:IMAS.SimpleSurface};
+                fixed_coils::Vector{Int}=Int[]) where {T<:Real}
     Nr, Nz = length(Rs), length(Zs)
     @assert size(Ψ) == (Nr, Nz)
     hr = Base.step(Rs)
@@ -163,10 +189,10 @@ function Canvas(Rs::AbstractRange{T},
     Vp  = zeros(length(surfaces))
     gm1 = zeros(length(surfaces))
     gm9 = zeros(length(surfaces))
-    r_cache, z_cache = IMASutils.contour_cache(Ψ)
+    r_cache, z_cache = IMASutils.contour_cache(Ψ; aggression_level=3)
     return Canvas(Rs, Zs, Ψ, Ip, coils, Rw, Zw, zt, zt, zt, zt, Ψpl, Ψvac, Gvac, Gbnd, U, Jt, Ψitp,
                   SVector{2,T}[], (0.0, 0.0), (0.0, 0.0), is_inside, is_in_wall, Rb_target, Zb_target,
-                  vs_circuit, rs_circuit, Ψ_at_coils, tmp_Ncoils, mutuals, mutuals_LU, a, b, c, MST, u,
+                  vs_circuit, rs_circuit, Ψ_at_coils, tmp_Ncoils, fixed_coils, mutuals, mutuals_LU, a, b, c, MST, u,
                   A, B, M, LU, S, tmp_Ψ, surfaces, Vp, gm1, gm9, r_cache, z_cache)
 end
 

src/flux_surfaces.jl

@@ -122,7 +122,7 @@ function compute_FSAs!(canvas::Canvas; update_surfaces=false)
     update_surfaces && trace_surfaces!(canvas)
     surfaces, Vp, gm1, gm9 = canvas._surfaces, canvas._Vp, canvas._gm1, canvas._gm9
 
-    sign_dpsi = sign(surfaces[end].psi - surfaces[1].psi)
+    sign_dpsi = (canvas.Ip >= 0.0) ? 1 : -1
 
     for (k, surface) in enumerate(surfaces)
 

src/workflow.jl

@@ -5,7 +5,7 @@ function solve!(canvas::Canvas, profile::CurrentProfile, Nout::Int, Nin::Int;
                 relax::Real=0.5,
                 tolerance::Real=0.0,
                 control::Union{Nothing, Symbol}=:shape,
-                fixed_coils::AbstractVector{Int}=Int[],
+                fixed_coils::AbstractVector{Int}=canvas._fixed_coils,
                 initialize_current=true,
                 initialize_mutuals=(control === :eddy))
 
@@ -47,7 +47,7 @@ function solve!(canvas::Canvas, profile::CurrentProfile, Nout::Int, Nin::Int;
                 @. Ψpl = (1.0 - relax) * Ψp0 + relax * Ψpl
             end
             update_bounds!(canvas; update_Ψitp=true)
-            Jtor!(canvas, profile; update_surfaces=false)
+            Jtor!(canvas, profile; update_surfaces=(j==1 && i==1))
             error_inner = abs((canvas.Ψaxis - Ψai) / (relax * Ψai))
             if sum(debug) == 2
                 println("\tInner $(i):\t", canvas.Ψaxis, "\t", canvas.Ψbnd - canvas.Ψaxis, "\t", error_inner)