Update to master branch

src/constants.jl

@@ -181,26 +181,27 @@ function Base.findall(identifier_name::Symbol, ids::IDSvector)
 end
 
 """
- Base.resize!(@nospecialize(ids::IDSvector{T}), identifier_name::Symbol; allow_multiple_matches=false)::T where {T<:IDSvectorElement}
+ Base.resize!(@nospecialize(ids::IDSvector{T}), identifier_name::Symbol, conditions::Pair{String}...; wipe::Bool=true, allow_multiple_matches::Bool=false)::T where {T<:IDSvectorElement}
 
-Resize ids if `identifier_name` is not found based on `ids.identifier.index` of `index_2_name(ids)`
+Resize ids if `identifier_name` is not found based on `ids.identifier.index` of `index_2_name(ids)` and a set of conditions are not met.
 
-If an entry matching the condition is found, then the content of the matching IDS is emptied, and the IDS is populated with the conditions.
+If wipe=true and an entry matching the condition is found, then the content of the matching IDS is emptied.
+
+Either way, the IDS is populated with the conditions.
 
 NOTE: `allow_multiple_matches` will delete all entries matching the conditions.
 
-Returns selected IDS(s)
+Returns the selected IDS
 """
-function Base.resize!(@nospecialize(ids::IDSvector{T}), identifier_name::Symbol, conditions::Pair{String}...; allow_multiple_matches=false)::T where {T<:IDSvectorElement}
+function Base.resize!(@nospecialize(ids::IDSvector{T}), identifier_name::Symbol, conditions::Pair{String}...; wipe::Bool=true, allow_multiple_matches::Bool=false)::T where {T<:IDSvectorElement}
  i = get(name_2_index(ids), identifier_name, nothing)
  if i === nothing
  error("`$(repr(identifier_name))` is not a known identifier for dd.$(fs2u(eltype(ids)))")
  elseif :identifier in fieldnames(eltype(ids))
- return resize!(ids, "identifier.index" => i, conditions...; allow_multiple_matches)
+ return resize!(ids, "identifier.index" => i, conditions...; wipe, allow_multiple_matches)
  else
- return resize!(ids, "index" => i, conditions...; allow_multiple_matches)
+ return resize!(ids, "index" => i, conditions...; wipe, allow_multiple_matches)
  end
- return eltype(haystack)[haystack[index] for index in indexes]
 end
 
 """

src/expressions/dynamic.jl

@@ -116,9 +116,6 @@ dyexp["core_profiles.profiles_1d[:].j_tor"] =
  Jpar_2_Jtor(rho_tor_norm, profiles_1d.j_total, true, eqt)
  end
 
-dyexp["core_profiles.profiles_1d[:].time"] =
- (; core_profiles, profiles_1d_index, _...) -> core_profiles.time[profiles_1d_index]
-
 # core_profiles.vacuum_toroidal_field #
 
 dyexp["core_profiles.vacuum_toroidal_field.b0"] =
@@ -144,9 +141,6 @@ dyexp["core_profiles.global_quantities.beta_tor_norm"] =
 #= ============ =#
 # core_transport #
 #= ============ =#
-dyexp["core_transport.model[:].profiles_1d[:].time"] =
- (; core_transport, profiles_1d_index, _...) -> core_transport.time[profiles_1d_index]
-
 dyexp["core_transport.vacuum_toroidal_field.b0"] =
  (time; dd, _...) -> vacuum_r0_b0_time(dd, time)[2]
 
@@ -268,9 +262,6 @@ dyexp["equilibrium.time_slice[:].profiles_1d.darea_dpsi"] =
  (psi; profiles_1d, _...) -> gradient(psi, profiles_1d.area)
 
 
-dyexp["equilibrium.time_slice[:].time"] =
- (; equilibrium, time_slice_index, _...) -> equilibrium.time[time_slice_index]
-
 dyexp["equilibrium.time_slice[:].profiles_1d.psi_norm"] =
  (psi; _...) -> norm01(psi)
 
@@ -394,13 +385,6 @@ dyexp["core_sources.source[:].profiles_1d[:].ion[:].particles"] =
  end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].time"] =
- (; core_sources, profiles_1d_index, _...) -> core_sources.time[profiles_1d_index]
-
-dyexp["core_sources.source[:].global_quantities[:].time"] =
- (; core_sources, global_quantities_index, _...) -> core_sources.time[global_quantities_index]
-
-
 dyexp["core_sources.vacuum_toroidal_field.b0"] =
  (time; dd, _...) -> vacuum_r0_b0_time(dd, time)[2]
 

src/math.jl

@@ -19,7 +19,7 @@ Turn a vector into a range (if possible)
 """
 function to_range(vector::AbstractVector{<:Real})
  tmp = diff(vector)
- if !(1 - sum(abs.(tmp .- tmp[1])) / length(vector) ≈ 1.0)
+ if !(1 - sum(abs, tmp .- tmp[1]) / length(vector) ≈ 1.0)
  error("to_range requires vector data to be equally spaced")
  end
  return range(vector[1], vector[end], length=length(vector))
@@ -354,7 +354,7 @@ end
 
 Calculate the angle between three points
 """
-function calculate_angle(p1::T, p2::T, p3::T) where {T<:AbstractVector{<:Real}}
+function calculate_angle(p1::T, p2::T, p3::T) where {T}
  v1 = [p2[1] - p1[1], p2[2] - p1[2]]
  v2 = [p3[1] - p2[1], p3[2] - p2[2]]
  dot_product = dot(v1, v2)
@@ -651,7 +651,22 @@ function angle_between_two_vectors(
  return acos((v1_x * v2_x + v1_y * v2_y) / (sqrt(v1_x^2 + v1_y^2) * sqrt(v2_x^2 + v2_y^2)))
 end
 
-function unique_indices(arr)
- uniq_elements = unique(arr)
- return [findfirst(==(elem), arr) for elem in uniq_elements]
-end
+"""
+ unique_indices(vec::AbstractVector)::Vector{Int}
+
+Return the indices of the first occurrence of each unique element in the input vector `vec`
+"""
+function unique_indices(vec::AbstractVector)::Vector{Int}
+ uniq_elements = unique(vec)
+ return [findfirst(==(elem), vec) for elem in uniq_elements]
+end
+
+"""
+ getindex_circular(vec::AbstractVector{T}, idx::Int)::T where {T}
+
+Return the element of the vector `vec` at the position `idx`. 
+If `idx` is beyond the length of `vec` or less than 1, it wraps around in a circular manner.
+"""
+function getindex_circular(vec::AbstractVector{T}, idx::Int)::T where {T}
+ return vec[(idx-1)%length(vec)+1]
+end

src/physics/build.jl

@@ -19,31 +19,15 @@ function get_build(bd::IMAS.build; kw...)
 end
 
 """
- get_build(
+ get_build_layers(
  layers::IMAS.IDSvector{<:IMAS.build__layer};
  type::Union{Nothing,BuildLayerType}=nothing,
  name::Union{Nothing,String}=nothing,
  identifier::Union{Nothing,Integer}=nothing,
- fs::Union{Nothing,BuildLayerSide,AbstractVector{BuildLayerSide}}=nothing,
- return_only_one::Bool=true,
- return_index::Bool=false,
- raise_error_on_missing::Bool=true)
+ fs::Union{Nothing,BuildLayerSide,AbstractVector{BuildLayerSide}}=nothing)
 
 Select layer(s) in build based on a series of selection criteria
-With `raise_error_on_missing=false` will returns `missing` if layer is missing
 """
-function get_build(
- layers::IMAS.IDSvector{<:IMAS.build__layer};
- type::Union{Nothing,BuildLayerType}=nothing,
- name::Union{Nothing,String}=nothing,
- identifier::Union{Nothing,Integer}=nothing,
- fs::Union{Nothing,BuildLayerSide,AbstractVector{BuildLayerSide}}=nothing,
- return_only_one::Bool=true,
- return_index::Bool=false,
- raise_error_on_missing::Bool=true)
- error("`IMAS.get_build()` is obsolete. Use `get_build_layer()` and `get_build_index()` or `get_build_layers()` and `get_build_indexes()` instead")
-end
-
 function get_build_layers(
  layers::IMAS.IDSvector{<:IMAS.build__layer};
  type::Union{Nothing,BuildLayerType}=nothing,
@@ -64,6 +48,16 @@ function get_build_layers(
  return valid_layers
 end
 
+"""
+ get_build_indexes(
+ layers::IMAS.IDSvector{<:IMAS.build__layer};
+ type::Union{Nothing,BuildLayerType}=nothing,
+ name::Union{Nothing,String}=nothing,
+ identifier::Union{Nothing,Integer}=nothing,
+ fs::Union{Nothing,BuildLayerSide,AbstractVector{BuildLayerSide}}=nothing)
+
+Returns indexes of layer(s) in build based on a series of selection criteria
+"""
 function get_build_indexes(
  layers::IMAS.IDSvector{<:IMAS.build__layer};
  type::Union{Nothing,BuildLayerType}=nothing,
@@ -84,6 +78,18 @@ function get_build_indexes(
  return valid_layers_indexes
 end
 
+"""
+ get_build_layer(
+ layers::IMAS.IDSvector{<:IMAS.build__layer};
+ type::Union{Nothing,BuildLayerType}=nothing,
+ name::Union{Nothing,String}=nothing,
+ identifier::Union{Nothing,Integer}=nothing,
+ fs::Union{Nothing,BuildLayerSide,AbstractVector{BuildLayerSide}}=nothing)
+
+Select layer in build based on a series of selection criteria
+
+It raises an error if none or more than one layer matches.
+"""
 function get_build_layer(
  layers::IMAS.IDSvector{<:IMAS.build__layer};
  type::Union{Nothing,BuildLayerType}=nothing,
@@ -102,6 +108,18 @@ function get_build_layer(
  return valid_layers[1]
 end
 
+"""
+ get_build_index(
+ layers::IMAS.IDSvector{<:IMAS.build__layer};
+ type::Union{Nothing,BuildLayerType}=nothing,
+ name::Union{Nothing,String}=nothing,
+ identifier::Union{Nothing,Integer}=nothing,
+ fs::Union{Nothing,BuildLayerSide,AbstractVector{BuildLayerSide}}=nothing)
+
+Returns index of layer in build based on a series of selection criteria
+
+It raises an error if none or more than one layer matches.
+"""
 function get_build_index(
  layers::IMAS.IDSvector{<:IMAS.build__layer};
  type::Union{Nothing,BuildLayerType}=nothing,
@@ -184,7 +202,7 @@ function structures_mask(bd::IMAS.build; ngrid::Int=257, border_fraction::Real=0
  end
 end
 
-function func_nested_layers(layer::IMAS.build__layer, func::Function)
+function func_nested_layers(layer::IMAS.build__layer{D}, func::Function)::D where {D<:Real}
  i = index(layer)
  layers = parent(layer)
  # _in_ layers or plasma

src/physics/equilibrium.jl

@@ -142,39 +142,39 @@ end
 """
  vacuum_r0_b0_time(dd::IMAS.dd) 
 
-Returns R0 as well as B0, and time arrays
+Returns `R0` as well as the `B0` and `time` arrays.
 
-This function solves the issue that in IMAS the information about R0 and B0 is replicated across different IDSs
+This function solves the issue that in IMAS the information about R0 and B0 is replicated across different IDSs.
 """
-function vacuum_r0_b0_time(dd::IMAS.dd)
+function vacuum_r0_b0_time(dd::IMAS.dd{T}) where {T<:Real}
  source = Set{Symbol}()
 
  # R0
  if hasfield(typeof(dd), :tf) && isfilled(dd.tf, :r0)
- R0 = dd.tf.r0
+ R0 = dd.tf.r0::T
  push!(source, :tf)
  else
  for (name, ids) in dd
  if hasfield(typeof(ids), :vacuum_toroidal_field)
  if isfilled(ids.vacuum_toroidal_field, :r0)
- R0 = ids.vacuum_toroidal_field.r0
+ R0 = ids.vacuum_toroidal_field.r0::T
  push!(source, name)
  break
  end
  end
  end
  end
 
- # B0 and time
+ # B0 and time vectors
  # from: tf
  if hasfield(typeof(dd), :tf) && isfilled(dd.tf.b_field_tor_vacuum_r, :data)
- B0 = dd.tf.b_field_tor_vacuum_r.data / R0
+ B0 = dd.tf.b_field_tor_vacuum_r.data::Vector{T} / R0
  time = dd.tf.b_field_tor_vacuum_r.data
  push!(source, :tf)
 
  # from: pulse_schedule
  elseif isfilled(dd.pulse_schedule.tf.b_field_tor_vacuum_r.reference, :data)
- B0 = dd.pulse_schedule.tf.b_field_tor_vacuum_r.reference.data / R0
+ B0 = dd.pulse_schedule.tf.b_field_tor_vacuum_r.reference.data::T / R0
  time = dd.pulse_schedule.tf.b_field_tor_vacuum_r.reference.time
  push!(source, :pulse_schedule)
 
@@ -185,8 +185,8 @@ function vacuum_r0_b0_time(dd::IMAS.dd)
  for (name, ids) in dd
  if hasfield(typeof(ids), :vacuum_toroidal_field)
  if isfilled(ids.vacuum_toroidal_field, :b0) && isfilled(ids, :time)
- B0_ = ids.vacuum_toroidal_field.b0
- time_ = ids.time
+ B0_ = ids.vacuum_toroidal_field.b0::Vector{T}
+ time_ = ids.time::Vector{Float64}
  append!(time, time_)
  append!(B0, B0_)
  reps = length(time_) - length(B0_)
@@ -199,6 +199,7 @@ function vacuum_r0_b0_time(dd::IMAS.dd)
  B0 = B0[index]
  time = time[index]
  end
+
  return R0, B0, time, source
 end
 

src/physics/nuclear.jl

@@ -223,7 +223,7 @@ function D_T_to_He4_source!(cs::IMAS.core_sources, cp::IMAS.core_profiles)
  reactivity = D_T_to_He4_reactions(cp1d; polarized_fuel_fraction)
  ion_to_electron_fraction = sivukhin_fraction(cp1d, eV1, 4.0)
  energy = reactivity .* eV1 * constants.e # J/m^3/s = W/m^3
- source = resize!(cs.source, :fusion, "identifier.name" => name)
+ source = resize!(cs.source, :fusion, "identifier.name" => name; wipe=false, allow_multiple_matches=true)
  new_source(
  source,
  source.identifier.index,
@@ -254,7 +254,7 @@ function D_D_to_He3_source!(cs::IMAS.core_sources, cp::IMAS.core_profiles)
  reactivity = D_D_to_He3_reactions(cp1d)
  ion_to_electron_fraction = sivukhin_fraction(cp1d, eV1, 3.0)
  energy = reactivity .* eV1 * constants.e # J/m^3/s = W/m^3
- source = resize!(cs.source, :fusion, "identifier.name" => name)
+ source = resize!(cs.source, :fusion, "identifier.name" => name; wipe=false, allow_multiple_matches=true)
  new_source(
  source,
  source.identifier.index,
@@ -284,7 +284,7 @@ function D_D_to_T_source!(cs::IMAS.core_sources, cp::IMAS.core_profiles)
  reactivity = D_D_to_T_reactions(cp1d)
  ion_to_electron_fraction = sivukhin_fraction(cp1d, eV1, 3.0)
  energy = reactivity .* eV1 * constants.e # J/m^3/s = W/m^3
- source = resize!(cs.source, :fusion, "identifier.name" => name)
+ source = resize!(cs.source, :fusion, "identifier.name" => name; wipe=false, allow_multiple_matches=true)
  new_source(
  source,
  source.identifier.index,
@@ -302,7 +302,7 @@ function D_D_to_T_source!(cs::IMAS.core_sources, cp::IMAS.core_profiles)
  eV2 = 3.0225e6
  ion_to_electron_fraction = sivukhin_fraction(cp1d, eV2, 1.0)
  energy = reactivity .* eV2 * constants.e # J/m^3/s = W/m^3
- source = resize!(cs.source, :fusion, "identifier.name" => name)
+ source = resize!(cs.source, :fusion, "identifier.name" => name; wipe=false, allow_multiple_matches=true)
  new_source(
  source,
  source.identifier.index,

src/physics/profiles.jl

@@ -189,9 +189,9 @@ function tau_e_thermal(cp1d::IMAS.core_profiles__profiles_1d, cs::IMAS.core_sour
 end
 
 
-function tau_e_h98(dd::IMAS.dd; time::Float64=dd.global_time)
- eqt = dd.equilibrium.time_slice[time]
- cp1d = dd.core_profiles.profiles_1d[time]
+function tau_e_h98(dd::IMAS.dd; time0::Float64=dd.global_time)
+ eqt = dd.equilibrium.time_slice[time0]
+ cp1d = dd.core_profiles.profiles_1d[time0]
  cs = dd.core_sources
  return tau_e_h98(eqt, cp1d, cs)
 end
@@ -223,9 +223,9 @@ function tau_e_h98(eqt::IMAS.equilibrium__time_slice, cp1d::IMAS.core_profiles__
  return tau98
 end
 
-function tau_e_ds03(dd::IMAS.dd; time::Float64=dd.global_time)
- eqt = dd.equilibrium.time_slice[time]
- cp1d = dd.core_profiles.profiles_1d[time]
+function tau_e_ds03(dd::IMAS.dd; time0::Float64=dd.global_time)
+ eqt = dd.equilibrium.time_slice[time0]
+ cp1d = dd.core_profiles.profiles_1d[time0]
  cs = dd.core_sources
  return tau_e_ds03(eqt, cp1d, cs)
 end

src/physics/radiation.jl

@@ -27,7 +27,7 @@ function bremsstrahlung_source!(dd::IMAS.dd)
 
  powerDensityBrem = -1.690e-38 .* ne .^ 2 .* cp1d.zeff .* sqrt.(Te)
 
- source = resize!(dd.core_sources.source, :bremsstrahlung)
+ source = resize!(dd.core_sources.source, :bremsstrahlung; wipe=false, allow_multiple_matches=true)
  new_source(source, source.identifier.index, "brem", cp1d.grid.rho_tor_norm, cp1d.grid.volume, cp1d.grid.area; electrons_energy=powerDensityBrem)
  return dd
 end
@@ -79,7 +79,7 @@ function synchrotron_source!(dd::IMAS.dd; wall_reflection_coefficient=0.8)
  # Synchrotron radiation
  powerDensitySync = rad_sync.(ϵ, a, B0, ne, Te; wall_reflection_coefficient)
 
- source = resize!(dd.core_sources.source, :synchrotron_radiation)
+ source = resize!(dd.core_sources.source, :synchrotron_radiation; wipe=false, allow_multiple_matches=true)
  new_source(source, source.identifier.index, "synch", cp1d.grid.rho_tor_norm, cp1d.grid.volume, cp1d.grid.area; electrons_energy=powerDensitySync)
  return source
 end
@@ -102,7 +102,7 @@ function line_radiation_source!(dd::IMAS.dd)
  linerad .+= rad_ion_adas(Te, ne, ni, zi, namei)
  end
 
- source = resize!(dd.core_sources.source, :line_radiation)
+ source = resize!(dd.core_sources.source, :line_radiation; wipe=false, allow_multiple_matches=true)
  new_source(source, source.identifier.index, "line", cp1d.grid.rho_tor_norm, cp1d.grid.volume, cp1d.grid.area; electrons_energy=linerad)
  return source
 end