use MXH δ ζ ovality twist

ProjectTorreyPines · Oct 5, 2024 · d432496 · d432496
1 parent 9170836
commit d432496
Show file tree

Hide file tree

Showing 4 changed files with 26 additions and 21 deletions.
diff --git a/src/actors/build/cx_actor.jl b/src/actors/build/cx_actor.jl
@@ -104,12 +104,11 @@ function segmented_wall(eq_r::AbstractVector{T}, eq_z::AbstractVector{T}, gap::T
 
         # automatic symmetry detection
         if symmetric === nothing
-            symmetric = sum(abs.(mxh.c)) / length(mxh.c) < 1E-3
+            symmetric = (abs(mxh.c0) .+ sum(abs.(mxh.c))) / (length(mxh.c) + 1) < 1E-3
         end
 
         # negative δ
-        δ = sin(mxh.s[1])
-        if δ < 0.0
+        if mxh.δ < 0.0
             Θ = LinRange(π / 4, 2 * π - π / 4, 100) # overshoot
         else
             Θ = LinRange(-π * 3 / 4, π * 3 / 4, 100) # overshoot
@@ -121,7 +120,7 @@ function segmented_wall(eq_r::AbstractVector{T}, eq_z::AbstractVector{T}, gap::T
         Z = (Z .- mxh.Z0) .* 1.1 .+ mxh.Z0
 
         # correct hfs to have a flat center stack wall
-        if δ < 0.0
+        if mxh.δ < 0.0
             R += IMAS.interp1d([Θ[1], Θ[argmax(Z)], Θ[argmin(Z)], Θ[end]], [maximum(eq_r) - R[1], 0.0, 0.0, maximum(eq_r) - R[end]]).(Θ)
         else
             R += IMAS.interp1d([Θ[1], Θ[argmax(Z)], Θ[argmin(Z)], Θ[end]], [minimum(eq_r) - R[1], 0.0, 0.0, minimum(eq_r) - R[end]]).(Θ)
@@ -327,7 +326,13 @@ function wall_from_eq!(
     return wall
 end
 
-function divertor_regions!(bd::IMAS.build{T}, eqt::IMAS.equilibrium__time_slice{T}, divertors::IMAS.divertors{T}, wall_r::AbstractVector{T}, wall_z::AbstractVector{T}) where {T<:Real}
+function divertor_regions!(
+    bd::IMAS.build{T},
+    eqt::IMAS.equilibrium__time_slice{T},
+    divertors::IMAS.divertors{T},
+    wall_r::AbstractVector{T},
+    wall_z::AbstractVector{T}
+) where {T<:Real}
     RA = eqt.global_quantities.magnetic_axis.r
     ZA = eqt.global_quantities.magnetic_axis.z
 
@@ -424,11 +429,11 @@ function divertor_regions!(bd::IMAS.build{T}, eqt::IMAS.equilibrium__time_slice{
         backwall_domain_poly = try
             LibGEOS.intersection(backwall_poly, domain_poly)
         catch e
-            display(plot(backwall_poly;aspect_ratio=:equal))
-            display(plot!(eqt;cx=true))
-            display(scatter!([RA,Rx],[ZA,Zx]))
-            display(plot!(xx,yy))
-            display(plot!(domain_poly;alpha=0.5))
+            display(plot(backwall_poly; aspect_ratio=:equal))
+            display(plot!(eqt; cx=true))
+            display(scatter!([RA, Rx], [ZA, Zx]))
+            display(plot!(xx, yy))
+            display(plot!(domain_poly; alpha=0.5))
             rethrow(e)
         end
         divertor_poly = LibGEOS.difference(backwall_domain_poly, plasma_poly)
@@ -634,7 +639,7 @@ function build_cx!(bd::IMAS.build, eqt::IMAS.equilibrium__time_slice, wall::IMAS
 
     # negative triangularity plasma
     mxh = IMAS.MXH(eqt.boundary.outline.r, eqt.boundary.outline.z, 1)
-    is_negative_D = sin(mxh.s[1]) < 0.0
+    is_negative_D = mxh.δ < 0.0
 
     #plot()
     plasma_to_tf = reverse(IMAS.get_build_indexes(bd.layer; fs=IMAS._hfs_))

diff --git a/src/ddinit/init_pulse_schedule.jl b/src/ddinit/init_pulse_schedule.jl
@@ -74,7 +74,7 @@ function init_pulse_schedule!(dd::IMAS.dd, ini::ParametersAllInits, act::Paramet
                 if ismissing(ini.rampup, :ends_at)
                     init_pulse_schedule_postion_control(pc, mxhb, time0)
                 else
-                    wr = wall_radii(mxhb.mxh.R0, mxhb.mxh.ϵ * mxhb.mxh.R0, ini.build.plasma_gap)
+                    wr = wall_radii(mxhb.mxh.R0, mxhb.mxh.minor_radius, ini.build.plasma_gap)
                     mxh_bore, mxh_lim2div = limited_to_diverted(0.75, mxhb, wr.r_hfs, wr.r_lfs, ini.rampup.side)
                     if time0 <= 0.0
                         init_pulse_schedule_postion_control(pc, mxh_bore, time0)
@@ -284,15 +284,15 @@ function init_pulse_schedule_postion_control(pc::IMAS.pulse_schedule__position_c
     IMAS.reorder_flux_surface!(pr, pz, argmax(pz))
 
     # scalars
-    IMAS.set_time_array(pc.minor_radius, :reference, time0, mxh.ϵ * mxh.R0)
+    IMAS.set_time_array(pc.minor_radius, :reference, time0, mxh.minor_radius)
     IMAS.set_time_array(pc.geometric_axis.r, :reference, time0, mxh.R0)
     IMAS.set_time_array(pc.geometric_axis.z, :reference, time0, mxh.Z0)
-    IMAS.set_time_array(pc.elongation, :reference, time0, mxh.κ)
-    IMAS.set_time_array(pc.tilt, :reference, time0, mxh.c0)
-    IMAS.set_time_array(pc.triangularity, :reference, time0, sin(mxh.s[1]))
-    IMAS.set_time_array(pc.squareness, :reference, time0, -mxh.s[2])
-    IMAS.set_time_array(pc.ovality, :reference, time0, mxh.c[1])
-    IMAS.set_time_array(pc.twist, :reference, time0, mxh.c[2])
+    IMAS.set_time_array(pc.elongation, :reference, time0, mxh.elongation)
+    IMAS.set_time_array(pc.tilt, :reference, time0, mxh.tilt)
+    IMAS.set_time_array(pc.triangularity, :reference, time0, mxh.triangularity)
+    IMAS.set_time_array(pc.squareness, :reference, time0, mxh.squareness)
+    IMAS.set_time_array(pc.ovality, :reference, time0, mxh.ovality)
+    IMAS.set_time_array(pc.twist, :reference, time0, mxh.twist)
 
     # boundary
     resize!(pc.boundary_outline, length(pr); wipe=false)

diff --git a/src/parameters/parameters_inits.jl b/src/parameters/parameters_inits.jl
@@ -600,7 +600,7 @@ Plots ini time dependent time traces including plasma boundary
 
         # plot equilibrium including x-points
         mxhb = MXHboundary(ini)
-        wr = wall_radii(mxhb.mxh.R0, mxhb.mxh.ϵ * mxhb.mxh.R0, ini.build.plasma_gap)
+        wr = wall_radii(mxhb.mxh.R0, mxhb.mxh.minor_radius, ini.build.plasma_gap)
         @series begin
             label := ""
             seriestype := :vline

diff --git a/src/physics.jl b/src/physics.jl
@@ -913,7 +913,7 @@ function fitMXHboundary(mxh::IMAS.MXH; upper_x_point::Bool, lower_x_point::Bool,
         if symmetric
             mxhb0.mxh.c0 = 0.0
             mxhb0.mxh.s = params[L+1:end]
-            mxhb0.mxh.c = mxhb0.mxh.s .* 0.0
+            mxhb0.mxh.c = zero(mxhb0.mxh.s)
         else
             L += 1
             mxhb0.mxh.c0 = params[L]