SOLPS2ctrl.jl

find_files_in_allowed_folders(
    input_dirs::String...;
    eqdsk_file::String,
    recursive::Bool=true,
)

Searches a list of allowed folders for a set of filenames that will provide information about the SOLPS case. Returns a list of filenames with complete paths.

Example:

SOLPS2ctrl.find_files_in_allowed_folders(
    "<your samples folder>/D3D_Ma_184833_03600";
    eqdsk_file="g184833.03600",
)

geqdsk_to_imas!(
    eqdsk_file::String,
    dd::IMAS.dd;
    set_time::Union{Nothing, Float64}=nothing,
    time_index::Int=1,
)

Transfers the equilibrium reconstruction from an EFIT-style gEQDSK file into the IMAS DD structure.

preparation(
    eqdsk_file::String,
    dirs::String...;
    core_method::String="simple",
    filename::String="sd_input_data",
    output_format::String="json",
    eqdsk_set_time::Union{Nothing, Float64}=nothing,
    eq_time_index::Int=1,
    allow_boundary_flux_correction::Bool=false,
)::IMAS.dd

Gathers SOLPS and EFIT files and loads them into IMAS structure. Extrapolates profiles as needed to get a complete picture.

function add_rho_to_equilibrium(dd:IMAS.dd)

Adds equilibrium rho profile to the DD

check_rho_1d(
    dd::IMAS.dd;
    time_slice::Int=1,
    throw_on_fail::Bool=false,
)::Bool

Checks to see if rho exists and is valid in the equilibrium 1d profiles

extrapolate_core(
    edge_rho::Vector{Float64},
    edge_quantity::Vector{Float64},
    rho_output::Vector{Float64},
)::Vector{Float64}

Function for assuming a core profile when given edge profile data.

Concept:

1. value and derivative should be continuous when joining real data
2. derivative at magnetic axis is known to be 0 when making profile vs. rho, by the definition of rho
3. derivative probably does something fancier between the pedestal and axis than just linear interpolation, so add an extra point in there
4. there's a joint between the steep pedestal and the shallow core that needs an extra knot to manage it properly
5. after making up a very simple gradient profile out of a few line segments, integrate it to get the profile of the quantity in question

fill_in_extrapolated_core_profile!(
    dd::IMAS.dd,
    quantity_name::String;
    method::String="simple",
    eq_time_idx::Int=1,
    eq_profiles_2d_idx::Int=1,
    value_field::Symbol=:values,
    grid_ggd_idx::Int=1,
    cell_subset_idx::Int=5,
)

This function accepts a DD that should be populated with equilibrium and edge_profiles as well as a request for a quantity to extrapolate into the core. It then maps edge_profiles data to rho, calls the function that performs the extrapolation (which is not a simple linear extrapolation but has some trickery to attempt to make a somewhat convincing profile shape), and writes the result to core_profiles. This involves a bunch of interpolations and stuff.

Input arguments:

• dd: an IMAS data dictionary
• quantity_name: the name of a quantity in edge_profiles.profiles_2d and core_profiles.profiles_1d, such as "electrons.density"
• method: Extrapolation method.
• eq_time_idx: index of the equilibrium time slice to use. For a typical SOLPS run, the SOLPS mesh will be based on the equilibrium reconstruction at a single time, so the DD associated with the SOLPS run only needs one equilibrium time slice to be loaded. However, one could combine the complete equilibrium time series with the SOLPS run and then have to specify which slice of the equilibrium corresponds to the SOLPS mesh.
• eq_profiles_2d_idx: index of the profiles_2D in equilibrium time_slice.
• value_field: Symbolic name of the values field in quantity.
• grid_ggd_idx: index of the grid_ggd to use. For a typical SOLPS run, the SOLPS grid is fixed, so this index defaults to 1. But in future, if a time varying grid is used, then this index will need to be specified.
• cell_subset_idx: index of the subset of cells to use for the extrapolation. The default is 5, which is the subset of all cells. If edge_profiles data is instead present for a different subset, for instance, -5, which are b2.5 cells only, then this index should be set to -5.

mesh_psi_spacing(
    dd::IMAS.dd;
    eq_time_idx::Int=1,
    eq_profiles_2d_idx::Int=1,
    grid_ggd_idx::Int=1,
    space_idx::Int=1,
    avoid_guard_cell::Bool=true,
    spacing_rule="mean",
)

Inspects the mesh to see how far apart faces are in psi_N. Requires that GGD and equilibrium are populated.

Input Arguments:

• dd: a data dictionary instance with required data loaded into it
• eq_time_idx: index of the equilibrium time slice to use. For a typical SOLPS run, the SOLPS mesh will be based on the equilibrium reconstruction at a single time, so the DD associated with the SOLPS run only needs one equilibrium time slice to be loaded. However, one could combine the complete equilibrium time series with the SOLPS run and then have to specify which slice of the equilibrium corresponds to the SOLPS mesh.
• eq_profiles_2d_idx: index of the profiles_2D in equilibrium time_slice.
• grid_ggd_idx: index of the grid_ggd to use. For a typical SOLPS run, the SOLPS grid is fixed, so this index defaults to 1. But in future, if a time varying grid is used, then this index will need to be specified.
• space_idx: index of the space to use. For a typical SOLPS run, there will be only one space so this index will mostly remain at 1.
• avoid_guard_cell: assume that the last cell is a guard cell so take end-2 and end-1 instead of end and end-1
• spacing_rule: "edge" or "mean" to make spacing of new cells (in psi_N) be the same as the spacing at the edge of the mesh, or the same as the average spacing

cached_mesh_extension!(
    dd::IMAS.dd,
    eqdsk_file::String,
    b2fgmtry::String;
    eq_time_idx::Int=1,
    eq_profiles_2d_idx::Int=1,
    grid_ggd_idx::Int=1,
    space_idx::Int=1,
    clear_cache::Bool=false,
)::String

Adds an extended mesh to a data dictionary, possibly from a cached result.

Input Arguments:

• dd: The data dictionary. It will be modified in place.
• eqdsk_file: the name of the EQDSK file that was used to get equilibrium data in the dd.
• b2fgmtry: the name of the SOLPS geometry file that was used to get GGD info in edge_profiles in the dd.
• eq_time_idx: Index of the time slice in equilibrium
• eq_profiles_2d_idx: Index of the 2D profile set in equilibrium (there is usually only one)
• grid_ggd_idx: Index of the grid_ggd set in edge_profiles
• space_idx: Index of the space
• clear_cache: delete any existing cache file (for use in testing)

offset_scale(
    val::Union{Float64, Vector{Float64}, Matrix{Float64}};
    offset::Union{Float64, Vector{Float64}}=0.0,
    factor::Union{Float64, Vector{Float64}}=1.0,
)::typeof(val)

Subtract an offset and multiply by a factor, the val to make it nominally in the range of -1 to 1 (not strictly) for easy identification of system.

(val .- offset) .* factor

unscale_unoffset(
    offset_scaled::Union{Float64, Vector{Float64}, Matrix{Float64}};
    offset::Union{Float64, Vector{Float64}}=0.0,
    factor::Union{Float64, Vector{Float64}}=1.0,
)::typeof(offset_scaled)

Undo previously applied offset and scaling.

offset_scaled ./ factor .+ offset

system_id(
    inp::Union{Vector{Float64}, Matrix{Float64}},
    out::Union{Vector{Float64}, Matrix{Float64}},
    tt::Vector{Float64},
    order::Int;
    input_cond::Union{Nothing, Function}=nothing,
    inp_offset::Float64=0.0, inp_factor::Float64=1.0,
    out_offset::Float64=0.0, out_factor::Float64=1.0,
    input_cond_kwargs::Dict{Symbol, T}=Dict{Symbol, Any}(),
    newpem_kwargs::Dict{Symbol, U}=Dict{Symbol, Any}(),
    verbose::Bool=false,
) where {T, U}

Perform system identification for a set on input data inp, output data out, and time series vector tt. If there are more than one inputs or outputs, provide them as Matrix with first dimension for ports (input or output) and second dimension for time. If input_cond is provided, it is applied to offseted and scaled input before performing system identification with keywords for this function provided in input_cond_kwargs.

This function uses ControlSystemIdentification.newpem to perform the system identification. Any additional keywords for this function should be passed as dictionary in newpem_kwargs. For advanced use, it is recommended to do system identification directly with newpem instead of using this function.

Returns a linear state space model of the order.

system_id_optimal_input_cond(
    inp::Union{Vector{Float64}, Matrix{Float64}},
    out::Union{Vector{Float64}, Matrix{Float64}},
    tt::Vector{Float64},
    order::Int,
    input_cond::Union{Nothing, Function},
    input_cond_args_guess::Dict{Symbol, T};
    inp_offset::Float64=0.0, inp_factor::Float64=1.0,
    out_offset::Float64=0.0, out_factor::Float64=1.0,
    input_cond_args_lower::Dict{Symbol, V}=Dict{Symbol, Any}(),
    input_cond_args_upper::Dict{Symbol, W}=Dict{Symbol, Any}(),
    newpem_kwargs::Dict{Symbol, U}=Dict{Symbol, Any}(),
    verbose::Bool=false,
) where {T, U, V, W}

Perform system identification for a set on input data inp, output data out, and time series vector tt. If there are more than one inputs or outputs, provide them as Matrix with first dimension for ports (input or output) and second dimension for time.

The input_cond is applied to offseted and scaled input before performing system identification. The input_cond_args_guess is used as initial keyword arguments that provide the parameters of the input_cond. These arguments are then used to find the best fit while iteratively performing system_id in each step.

This function uses ControlSystemIdentification.newpem to perform the system identification. Any additional keywords for this function should be passed as dictionary in newpem_kwargs.

This function uses LsqFit.curve_fit to fit the parameters of input conditions along with performing the system identification.

For advanced use, it is recommended to do system identification directly with newpem and optimize using your favorite fitting method instead of using this function.

Returns a linear state space model of the order and the keyword argument dictionary containing optimal parameters for input_cond function.

model_evolve(
    sys::Union{PredictionStateSpace, StateSpace},
    inp::Union{Vector{Float64}, Matrix{Float64}};
    input_cond::Union{Nothing, Function}=nothing,
    inp_offset::Float64=0.0, inp_factor::Float64=1.0,
    out_offset::Float64=0.0, out_factor::Float64=1.0,
    input_cond_kwargs::Dict{Symbol, T}=Dict{Symbol, Any}(),
)::Union{Vector{Float64}, Matrix{Float64}} where {T}

Evolve a state space model sys with the input steps inp. The input is offseted and scaled with inp_offset and inp_factor and the output is unscaled and unoffseted with out_offset and out_factor. If a function is provided as inp_cond it is applied to the input after scaling and offseting along with any keywords passed for it.

gas_unit_converter(
    value_in::Float64,
    units_in::String,
    units_out::String;
    species::String="H",
    temperature::Float64=293.15,
)

Converts gas flows between different units. Uses ideal gas law to convert between Pressure * volume type flows / quantities and count / current types of units. There is a version that accepts floats in and outputs floats, and another that deals in Unitful quantities.

gas_unit_converter(
    value_in::Unitful.Quantity,
    units_in::String,
    units_out::String;
    species::String="H",
    temperature=293.15 * Unitful.K,
)

Converts gas flows between different units. Uses ideal gas law to convert between Pressure * volume type flows / quantities and count / current types of units. This is the Unitful version.

Output will be unitful, but the units are not simplified automatically. You can perform operations such as

• (output |> Unitful.upreferred).val
• Unitful.uconvert(Unitful.whatever, output).val

to handle simplification or conversion of units.

Although this function pretends torr L s$^{-1}$ and Pa m$^3$ s$^{-1}$ are different, use of Unitful should cause them to behave the same way as long as you simplify or convert units at the end. This means that you can use other pressure*volume type gas units and call them torr L s$^{-1}$ and the script will deal with them up to having messy units in the output.