Merge pull request #48 from ProjectTorreyPines/fix_bcp_10-29

Update boundary_control_points functions

Author: orso82

src/control.jl

@@ -134,22 +134,44 @@ function IsoControlPoints(Rs::AbstractVector{T}, Zs::AbstractVector{T}) where {T
 end
 
 """
-    boundary_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999; Npts::Integer=99)
+    boundary_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999, ψbound::Real=0.0; Npts::Integer=99)
+Return a Vector of FluxControlPoint, each with target `ψbound`, at `Npts` equally distributed `fraction_inside` percent inside the the boundary of `EQfixed`
+"""
+function boundary_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999, ψbound::Real=0.0; Npts::Integer=99)
+    ψ0, ψb = MXHEquilibrium.psi_limits(EQfixed)
+    Sp = MXHEquilibrium.flux_surface(EQfixed, fraction_inside * (ψb - ψ0) + ψ0; n_interp=Npts)
+    ψtarget = fraction_inside * (ψb - ψ0) + ψ0 - ψb + ψbound
+    return [FluxControlPoint(Sp.r[k], Sp.z[k], ψtarget, 1.0 / Npts) for k in 1:length(Sp.r)-1]
+end
+
+"""
+    boundary_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999, ψbound::Real=0.0; Npts::Integer=99)
+Return a Vector of FluxControlPoint, each with target `ψbound`, at `Npts` equally distributed `fraction_inside` percent inside the the boundary of `shot`
+"""
+function boundary_control_points(shot::TEQUILA.Shot, fraction_inside::Float64=0.999, ψbound::Real=0.0; Npts::Integer=99)
+    bnd = MillerExtendedHarmonic.MXH(shot, fraction_inside)
+    θs = LinRange(0, 2π, Npts + 1)
+    ψtarget = ψbound + TEQUILA.psi_ρθ(shot, fraction_inside, 0.0)
+    return [FluxControlPoint(bnd(θ)..., ψtarget, 1.0 / Npts) for θ in θs[1:end-1]]
+end
+
+"""
+    boundary_iso_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999; Npts::Integer=99)
 
 Return a Vector of IsoControlPoints, at `Npts` equally distributed `fraction_inside` percent inside the the boundary of `EQfixed`
 """
-function boundary_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999; Npts::Integer=99)
+function boundary_iso_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999; Npts::Integer=99)
     ψ0, ψb = MXHEquilibrium.psi_limits(EQfixed)
     Sp = MXHEquilibrium.flux_surface(EQfixed, fraction_inside * (ψb - ψ0) + ψ0; n_interp=Npts)
     return VacuumFields.IsoControlPoints(Sp.r, Sp.z)
 end
 
 """
-    boundary_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999; Npts::Integer=99)
+    boundary_iso_control_points(EQfixed::MXHEquilibrium.AbstractEquilibrium, fraction_inside::Float64=0.999; Npts::Integer=99)
 
 Return a Vector of IsoControlPoints, at `Npts` equally distributed `fraction_inside` percent inside the the boundary of `shot`
 """
-function boundary_control_points(shot::TEQUILA.Shot, fraction_inside::Float64=0.999; Npts::Integer=99)
+function boundary_iso_control_points(shot::TEQUILA.Shot, fraction_inside::Float64=0.999; Npts::Integer=99)
     bnd = MillerExtendedHarmonic.MXH(shot, fraction_inside)
     θs = LinRange(0, 2π, Npts + 1)
     r = [bnd(θ)[1] for θ in θs[1:end-1]]

test/runtests.jl

@@ -130,7 +130,8 @@ end
         println("Testing for COCOS ", cc)
         Gfixed = efit(transform_cocos(gfixed, cc0, cc), cc)
         Gfree = efit(transform_cocos(gfree, cc0, cc), cc)
-        flux_cps = boundary_control_points(Gfree, 0.999)
+        _, ψbound = psi_limits(Gfree)
+        flux_cps = boundary_control_points(Gfree, 0.999, ψbound)
         saddle_cps = [SaddleControlPoint(Rx, Zx)]
         currents, _ = find_coil_currents!(coils, Gfixed; flux_cps, saddle_cps, λ_regularize=1e-14)
         if cc < 10
@@ -159,12 +160,15 @@ end
     Gfixed = efit(transform_cocos(gfixed, cc0, 11), 11)
     Gfree = efit(transform_cocos(gfree, cc0, 11), 11)
 
+    _, ψbound = psi_limits(Gfree)
+    Rs, Zs = gfree.r, gfree.z
+
     ix = argmin(gfree.zbbbs)
     Rx, Zx = gfree.rbbbs[ix], gfree.zbbbs[ix]
 
-    flux_cps = boundary_control_points(Gfree, 0.999)
+    flux_cps = boundary_control_points(Gfree, 0.999, ψbound)
     saddle_cps = [SaddleControlPoint(Rx, Zx)]
-    fixed2free(Gfixed, coils, gfree.r, gfree.z; flux_cps, saddle_cps);
+    fixed2free(Gfixed, coils, Rs, Zs; flux_cps, saddle_cps)
 end
 
 @testset "current_BtIp" begin
@@ -192,17 +196,23 @@ end
 
         gfree = MXHEquilibrium.readg(EQs_free[i]; set_time=0.0)
         Gfree = MXHEquilibrium.efit(gfree, cc)
+        _, ψbound = psi_limits(Gfree)
 
         #  Currents from EFIT
         p = plot(C[i, :]; linewidth=3, linecolor=:black)
 
         # Currents with ψbound=0
-        flux_cps = boundary_control_points(Gfixed, 0.999)
+        flux_cps = boundary_control_points(Gfixed, 0.999, 0.0)
         c0, _ = find_coil_currents!(coils, Gfixed; flux_cps, λ_regularize=1e-14)
 
         plot!(c0; linewidth=3, linecolor=:red)
 
+        # Currents with ψbound from EFIT
+        flux_cps = boundary_control_points(Gfixed, 0.999, ψbound)
+        cb, _ = find_coil_currents!(coils, Gfixed; flux_cps, λ_regularize=1e-14)
+        plot!(cb; linewidth=3, linecolor=:blue)
         display(p)
+
     end
 end