Skip to content

Apply averages to coil and magnetic data #71

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Apply averages to coil and magnetic data #71

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
24 changes: 11 additions & 13 deletions src/coils.jl
View file
Open in desktop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't think we should change the default behavior to be averaging. That's a specific use case. Often you just want the current at a specific time.

Original file line number Diff line number Diff line change
Expand Up @@ -105,18 +105,9 @@ elements(coil::Union{IMAScoil,IMASloop}) = Iterators.filter(!isempty, coil.eleme

# BCL 2/27/24
# N.B.: Not sure about sign with turns and such
function current_per_turn(coil::IMAScoil)
if ismissing(coil.current, :data)
return 0.0
end
return IMAS.@ddtime(coil.current.data)
end

function current_per_turn(loop::IMASloop)
if ismissing(loop, :current)
return 0.0
end
return IMAS.@ddtime(loop.current)
function current_per_turn(pf::Union{IMAScoil,IMASloop})
avg = IMAS.getavg(pf, IMAS.global_time(pf), 0.005)
return avg.data.val
end

@inline function set_current_per_turn!(coil::AbstractSingleCoil, current_per_turn::Real)
Expand Down Expand Up @@ -514,7 +505,14 @@ end

max_current(coil::AbstractCoil) = abs(current_per_turn(coil) * turns(coil))
max_current(mcoil::MultiCoil) = isempty(mcoil.coils) ? 0.0 : maximum(abs, current_per_turn(coil) * turns(coil) for coil in mcoil.coils)
max_current(icoil::IMAScoil) = ismissing(icoil.current, :data) ? 0.0 : (abs(IMAS.@ddtime(icoil.current.data) * maximum(turns(elm) for elm in icoil.element)))
function max_current(icoil::IMAScoil)
if ismissing(icoil.current, :data)
return 0.0
else
avg = IMAS.getavg(icoil, IMAS.global_time(icoil), 0.005)
return (abs(avg.data.val * maximum(turns(elm) for elm in icoil.element)))
end
end

@recipe function plot_coils(coils::AbstractVector{<:AbstractCoil}; color_by=:current, cname=:diverging)
if !isempty(coils)
Expand Down
31 changes: 18 additions & 13 deletions src/imas_control.jl
View file
Open in desktop
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Here I'm more fine with averaging for the flux loop since that's actual data that could be noisy. I suppose we should make sure that this doesn't have any issue with smooth, synthetic data or any need for causal averaging.

Original file line number Diff line number Diff line change
Expand Up @@ -32,19 +32,21 @@ function magnetic_control_points(dd::IMAS.dd{T};
magnetic_probe_weights::AbstractVector{<:Real}=T[]) where {T<:Real}

mags = dd.magnetics
pavgs = IMAS.getavg(Val(:b_field_pol_probe), dd, dd.global_time, 0.005)
pdata = pavgs.data

# Probe control points (Note: type and size ignored)
min_σ = minimum(IMAS.@ddtime(probe.field.data_σ) for probe in mags.b_field_pol_probe)
min_σ = minimum(probe.err for probe in pdata)
field_cps = VacuumFields.FieldControlPoint{T}[]
for (k, probe) in enumerate(mags.b_field_pol_probe)
!isempty(probe.field.validity) && probe.field.validity < 0 && continue
weight = isempty(magnetic_probe_weights) ? 1.0 : magnetic_probe_weights[k]
if !isempty(probe.field.data_σ)
IMAS.@ddtime(probe.field.data_σ) < eps() && continue
weight /= IMAS.@ddtime(probe.field.data_σ) / min_σ
pdata[k].err < eps() && continue
weight /= pdata[k].err / min_σ
end
#IMAS allows probes to mix toroidal and poloidal fields so that needs to be accounted for
bpol_field = isempty(probe.toroidal_angle) ? IMAS.@ddtime(probe.field.data) : IMAS.@ddtime(probe.field.data) * (1 - cos(probe.poloidal_angle) * sin(probe.toroidal_angle))
bpol_field = isempty(probe.toroidal_angle) ? pdata[k].val : pdata[k].val * (1 - cos(probe.poloidal_angle) * sin(probe.toroidal_angle))

push!(field_cps, VacuumFields.FieldControlPoint{T}(probe.position.r, probe.position.z, probe.poloidal_angle, bpol_field, weight))
end
Expand All @@ -53,19 +55,22 @@ function magnetic_control_points(dd::IMAS.dd{T};

iref = reference_flux_loop_index
loops = mags.flux_loop
flavgs = IMAS.getavg(Val(:flux_loop), dd, dd.global_time, 0.005)
fldata = flavgs.data

flux_cps = VacuumFields.FluxControlPoint{T}[]
loop_cps = VacuumFields.IsoControlPoint{T}[]
if iref > 0
# Fit a single reference flux value and the difference between the other flux_loops and this value (this matches how data is collected on DIII-D and fit with EFIT)
N = length(loops)
@assert iref <= N
# Convert from total flux uncertainty (IMAS convention) to differential
relative_σ(k) = sqrt(IMAS.@ddtime(loops[k].flux.data_σ)^2 - IMAS.@ddtime(loops[iref].flux.data_σ)^2)
relative_σ(k) = sqrt(fldata[k].err^2 - fldata[iref].err^2)
min_σ = 1 / eps()
for k in (iref+1):(iref+N-1) # exclude i_ref explicitly
ck = IMAS.index_circular(N, k)
if !isempty(loops[ck].flux.data_σ)
IMAS.@ddtime(loops[ck].flux.data_σ) < eps() && continue
fldata[ck].err < eps() && continue
min_σ = relative_σ(ck) < min_σ ? relative_σ(ck) : min_σ
end
end
Expand All @@ -75,7 +80,7 @@ function magnetic_control_points(dd::IMAS.dd{T};
!isempty(loops[ck].flux.validity) && loops[ck].flux.validity < 0 && continue
weight = isempty(flux_loop_weights) ? 1.0 : flux_loop_weights[ck]
if !isempty(loops[ck].flux.data_σ)
IMAS.@ddtime(loops[ck].flux.data_σ) < eps() && continue
fldata[ck].err < eps() && continue
weight /= relative_σ(ck) / min_σ
end

Expand All @@ -86,7 +91,7 @@ function magnetic_control_points(dd::IMAS.dd{T};
loops[ck].position[1].z,
loops[iref].position[1].r,
loops[iref].position[1].z,
IMAS.@ddtime(loops[ck].flux.data) - IMAS.@ddtime(loops[iref].flux.data),
fldata[ck].val - fldata[iref].val,
weight
)
)
Expand All @@ -96,19 +101,19 @@ function magnetic_control_points(dd::IMAS.dd{T};
@assert loops[iref].flux.validity >= 0
end

push!(flux_cps, VacuumFields.FluxControlPoint{T}(loops[iref].position[1].r, loops[iref].position[1].z, IMAS.@ddtime(loops[iref].flux.data), weight))
push!(flux_cps, VacuumFields.FluxControlPoint{T}(loops[iref].position[1].r, loops[iref].position[1].z, fldata[iref].val, weight))
else
# Fit each flux loop separately (this is how data is collected on NSTX and fit with EFIT)
min_σ = minimum(IMAS.@ddtime(loop.flux.data_σ) for loop in loops)
min_σ = minimum(loop.err for loop in fldata)
for (k, loop) in enumerate(loops)
!isempty(loop.flux.validity) && loop.flux.validity < 0 && continue
weight = isempty(flux_loop_weights) ? 1.0 : flux_loop_weights[k]
if !isempty(loop.flux.data_σ)
IMAS.@ddtime(loop.flux.data_σ) < eps() && continue
weight /= IMAS.@ddtime(loop.flux.data_σ) / min_σ
fldata[k].err < eps() && continue
weight /= fldata[k].err / min_σ
end

push!(flux_cps, VacuumFields.FluxControlPoint{T}(loop.position[1].r, loop.position[1].z, IMAS.@ddtime(loop.flux.data), weight))
push!(flux_cps, VacuumFields.FluxControlPoint{T}(loop.position[1].r, loop.position[1].z, fldata[k].val, weight))
end
end
return flux_cps, loop_cps, field_cps
Expand Down