diff --git a/Project.toml b/Project.toml index 980dc5f..a1d8b01 100644 --- a/Project.toml +++ b/Project.toml @@ -1,16 +1,21 @@ name = "SOLPS2ctrl" uuid = "a531d12f-ac8a-43e8-b6d9-bd121431dd49" authors = ["David Eldon "] -version = "2.1.0" +version = "2.2.0" [deps] ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" +ControlSystemIdentification = "3abffc1c-5106-53b7-b354-a47bfc086282" +ControlSystemsBase = "aaaaaaaa-a6ca-5380-bf3e-84a91bcd477e" +DelimitedFiles = "8bb1440f-4735-579b-a4ab-409b98df4dab" EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" IMAS = "13ead8c1-b7d1-41bb-a6d0-5b8b65ed587a" IMASggd = "b7b5e640-9b39-4803-84eb-376048795def" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +LsqFit = "2fda8390-95c7-5789-9bda-21331edee243" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" @@ -23,18 +28,23 @@ YAML = "ddb6d928-2868-570f-bddf-ab3f9cf99eb6" [compat] ArgParse = "1" -Contour = "0.6.3" -EFIT = "1.0.0" -IMAS = "1.2.2, 2, 3" -IMASggd = "3.1" -Interpolations = "0.15.1" -JSON = "0.21.4" -PhysicalConstants = "0.2.3" -PlotUtils = "1.4.1" -Plots = "1.40.3" -PolygonOps = "0.1.2" +Contour = "0" +ControlSystemIdentification = "2" +ControlSystemsBase = "1" +DelimitedFiles = "1" +EFIT = "1" +IMAS = "1, 2, 3, 4, 5" +IMASggd = "3.3" +Interpolations = "0" +JSON = "0" +LinearAlgebra = "1" +LsqFit = "0.15.1" +PhysicalConstants = "0" +PlotUtils = "1" +Plots = "1" +PolygonOps = "0" SOLPS2imas = "2.2" -Statistics = "1.9.0" -Unitful = "1.20.0" -YAML = "0.4.11" +Statistics = "1" +Unitful = "1" +YAML = "0" julia = "1.10" diff --git a/docs/src/index.md b/docs/src/index.md index 80f90e1..6f76813 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -71,6 +71,16 @@ mesh_psi_spacing cached_mesh_extension! ``` +## System identification and modeling + +```@docs +offset_scale +unscale_unoffset +system_id +system_id_optimal_input_cond +model_evolve +``` + ## Unit conversion utilities ```@docs diff --git a/sample/D3D_Lore_Time_Dep/.gitignore b/sample/D3D_Lore_Time_Dep/.gitignore new file mode 100644 index 0000000..5d433b1 --- /dev/null +++ b/sample/D3D_Lore_Time_Dep/.gitignore @@ -0,0 +1,2 @@ +/imas_with_ifo.json +/input_gas_atomsps.txt diff --git a/sample/D3D_Lore_Time_Dep/imas_with_ifo.json.dvc b/sample/D3D_Lore_Time_Dep/imas_with_ifo.json.dvc new file mode 100644 index 0000000..55ac2c9 --- /dev/null +++ b/sample/D3D_Lore_Time_Dep/imas_with_ifo.json.dvc @@ -0,0 +1,12 @@ +md5: d7a572fef4d1491734ca69cc9ad6f0fd +frozen: true +deps: +- path: D3D_Lore_Time_Dep/imas_with_ifo.json + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: 34463adc4d9b452712c39a52654dd8894af37820 +outs: +- md5: 8198d08d5dd62b797901546e58794106 + size: 249503911 + hash: md5 + path: imas_with_ifo.json diff --git a/sample/D3D_Lore_Time_Dep/input_gas_atomsps.txt.dvc b/sample/D3D_Lore_Time_Dep/input_gas_atomsps.txt.dvc new file mode 100644 index 0000000..b744005 --- /dev/null +++ b/sample/D3D_Lore_Time_Dep/input_gas_atomsps.txt.dvc @@ -0,0 +1,12 @@ +md5: 625ffa07ec8b07be8d65c190fc7e6493 +frozen: true +deps: +- path: D3D_Lore_Time_Dep/input_gas_atomsps.txt + repo: + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: 34463adc4d9b452712c39a52654dd8894af37820 +outs: +- md5: d731f64bcce6c9e656d11513503fd72f + size: 25490 + hash: md5 + path: input_gas_atomsps.txt diff --git a/src/SOLPS2ctrl.jl b/src/SOLPS2ctrl.jl index 1044aee..04cc4fa 100644 --- a/src/SOLPS2ctrl.jl +++ b/src/SOLPS2ctrl.jl @@ -10,6 +10,7 @@ export find_files_in_allowed_folders, geqdsk_to_imas!, preparation include("$(@__DIR__)/supersize_profile.jl") include("$(@__DIR__)/repair_eq.jl") include("$(@__DIR__)/unit_utils.jl") +include("$(@__DIR__)/system_id.jl") """ find_files_in_allowed_folders( @@ -198,8 +199,8 @@ function geqdsk_to_imas!( if nprim > 0 bsx = eqt.boundary_separatrix.x_point resize!(bsx, nprim) - xrprim = xrs[xseps.==1] - xzprim = xzs[xseps.==1] + xrprim = xrs[xseps .== 1] + xzprim = xzs[xseps .== 1] for i ∈ nprim bsx[i].r = xrprim[i] bsx[i].z = xzprim[i] @@ -209,8 +210,8 @@ function geqdsk_to_imas!( if nsec > 0 bssx = eqt.boundary_secondary_separatrix.x_point resize!(bssx, nsec) - xrsec = xrs[xseps.==2] - xzsec = xzs[xseps.==2] + xrsec = xrs[xseps .== 2] + xzsec = xzs[xseps .== 2] for i ∈ nsec bssx[i].r = xrsec[i] bssx[i].z = xzsec[i] @@ -311,7 +312,7 @@ function preparation( for core_profile ∈ core_profiles tags = split(core_profile, ".") parent = dd.edge_profiles.ggd[1] - for tag ∈ tags[1:end-1] + for tag ∈ tags[1:(end-1)] parent = getproperty(parent, Symbol(tag)) end qty = getproperty(parent, Symbol(tags[end]), core_profile) diff --git a/src/repair_eq.jl b/src/repair_eq.jl index 05af026..36fd9b3 100644 --- a/src/repair_eq.jl +++ b/src/repair_eq.jl @@ -169,7 +169,7 @@ function add_rho_to_equilibrium!(dd::IMAS.dd) z1i = Interpolations.linear_interpolation(r1, z1) z2i = Interpolations.linear_interpolation(r2, z2) rr = LinRange(rmin, rmax, 101) - rc = (rr[1:end-1] + rr[2:end]) / 2.0 + rc = (rr[1:(end-1)] + rr[2:end]) / 2.0 integral_part_ = [ log(rr[i+1] / rr[i]) * abs(z1i(rc[i]) - z2i(rc[i])) for i ∈ 1:length(rc) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 09162f8..45fee4b 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -5,7 +5,8 @@ Utilities for extrapolating profiles using Interpolations: Interpolations using IMASggd: IMASggd, get_grid_subset, add_subset_element!, get_subset_boundary, - project_prop_on_subset!, get_subset_centers, get_TPS_mats, subset_do + project_prop_on_subset!, get_subset_centers, get_TPS_mats, subset_do, + all__grid_subset_prop using PolygonOps: PolygonOps using JSON: JSON @@ -91,12 +92,14 @@ function extrapolate_core( q_extend .+= q_offset output_profile = Array{Float64}(undef, length(rho_output)) - output_profile[rho_output.=rf] = + output_profile[rho_output .< rf] = q_extend[rho_output .< rf] + output_profile[rho_output .>= rf] = Interpolations.linear_interpolation( edge_rho, edge_quantity, - ).(rho_output[rho_output.>=rf]) + ).( + rho_output[rho_output .>= rf], + ) return output_profile end @@ -107,8 +110,9 @@ end method::String="simple", eq_time_idx::Int=1, eq_profiles_2d_idx::Int=1, + value_field::Symbol=:values, grid_ggd_idx::Int=1, - space_idx::Int=1, + cell_subset_idx::Int=5, ) This function accepts a DD that should be populated with `equilibrium` and @@ -131,11 +135,10 @@ Input arguments: 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. - - `space_id`x: 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. - `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 @@ -147,13 +150,12 @@ function fill_in_extrapolated_core_profile!( method::String="simple", eq_time_idx::Int=1, eq_profiles_2d_idx::Int=1, + value_field::Symbol=:values, grid_ggd_idx::Int=1, - space_idx::Int=1, cell_subset_idx::Int=5, ) check_rho_1d(dd; time_slice=eq_time_idx, throw_on_fail=true) grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] - space = grid_ggd.space[space_idx] cell_subset = get_grid_subset(grid_ggd, cell_subset_idx) midplane_subset = get_grid_subset(grid_ggd, 11) @@ -183,20 +185,17 @@ function fill_in_extrapolated_core_profile!( resize!(dd.core_profiles.profiles_1d, nt) end TPS_mats = get_TPS_mats(grid_ggd, cell_subset_idx) - for it ∈ 1:nt - tags = split(quantity_name, ".") - quantity_str = dd.edge_profiles.ggd[it] - for tag ∈ tags - quantity_str = getproperty(quantity_str, Symbol(tag)) - end - midplane_cell_centers, quantity = project_prop_on_subset!( - quantity_str, - cell_subset, - midplane_subset, - space; - TPS_mats=TPS_mats, - ) + mcc_qty = project_prop_on_subset!( + dd.edge_profiles.ggd, + quantity_name, + cell_subset, + midplane_subset; + value_field, + TPS_mats) + + for it ∈ 1:nt + midplane_cell_centers, quantity = mcc_qty[it] # Now quantity is at the outboard midplane # Get the rho values to go with the midplane quantity values @@ -239,7 +238,9 @@ function fill_in_extrapolated_core_profile!( Interpolations.linear_interpolation( psi1_eq, rho1_eq, - ).(psi_for_quantity[in_bounds]) + ).( + psi_for_quantity[in_bounds], + ) # Make sure the output 1D rho grid exists; create it if needed if IMAS.ismissing(dd.core_profiles.profiles_1d[it].grid, :rho_tor_norm) @@ -260,7 +261,7 @@ function fill_in_extrapolated_core_profile!( end parent = dd.core_profiles.profiles_1d[it] tags = split(quantity_name, ".") - for tag ∈ tags[1:end-1] + for tag ∈ tags[1:(end-1)] parent = getproperty(parent, Symbol(tag)) end setproperty!(parent, Symbol(tags[end]), quantity_core) @@ -415,7 +416,7 @@ function mesh_psi_spacing( end else if avoid_guard_cell - dpsin = diff(psin[2:end-1]) + dpsin = diff(psin[2:(end-1)]) else dpsin = diff(psin) end diff --git a/src/system_id.jl b/src/system_id.jl new file mode 100644 index 0000000..c139ea2 --- /dev/null +++ b/src/system_id.jl @@ -0,0 +1,279 @@ +import ControlSystemsBase: lsim, StateSpace +import ControlSystemIdentification: iddata, newpem, PredictionStateSpace +import LinearAlgebra: I, inv +import LsqFit: curve_fit, coef + +export offset_scale, + unscale_unoffset, model_evolve, system_id, system_id_optimal_input_cond +""" + 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 +""" +function 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) + return (val .- offset) .* factor +end + +""" + 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 +""" +function 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) + return offset_scaled ./ factor .+ offset +end + +""" + 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](https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss/#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`. +""" +function 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} + @assert size(inp)[end] == size(out)[end] == length(tt) + inp_so = offset_scale(inp; offset=inp_offset, factor=inp_factor) + if !isnothing(input_cond) + inp_so = input_cond(inp_so; input_cond_kwargs...) + end + u = Matrix{Float64}[] + if isa(inp_so, Vector) + u = Matrix(inp_so') + else + u = inp_so + end + + out_so = offset_scale(out; offset=out_offset, factor=out_factor) + v = Matrix{Float64}[] + if isa(out_so, Vector) + v = Matrix(out_so') + else + v = out_so + end + + id_data = iddata(v, u, tt[2] - tt[1]) + + # Following suppression of stdout and stderr is because of this issue + # in newpem. + # https://github.com/baggepinnen/ControlSystemIdentification.jl/issues/184 + # Once this issue has been resolved, we should only suppress stdout + original_stdout = stdout + original_stderr = stderr + if !verbose + redirect_stdout(devnull) + redirect_stderr(devnull) + end + sys = newpem(id_data, order; newpem_kwargs...).sys + if !verbose + redirect_stdout(original_stdout) + redirect_stderr(original_stderr) + end + return sys +end + +""" + 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. +""" +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}(), +)::typeof(inp) where {T} + inp_so = offset_scale(inp; offset=inp_offset, factor=inp_factor) + if !isnothing(input_cond) + inp_so = input_cond(inp_so; input_cond_kwargs...) + end + u = Matrix{Float64}[] + if isa(inp_so, Vector) + u = Matrix(inp_so') + else + u = inp_so + end + + x0 = inv(Matrix{Float64}(I, size(sys.A)...) - sys.A) * sys.B * u[:, 1] + modeled_out_so, _, x0, _ = lsim(sys, u; x0) + modeled_out = unscale_unoffset(modeled_out_so; offset=out_offset, factor=out_factor) + if size(modeled_out)[1] == 1 + return modeled_out[1, :] + else + return modeled_out + end +end + +""" + 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`](@ref) in each step. + +This function uses [ControlSystemIdentification.newpem](https://baggepinnen.github.io/ControlSystemIdentification.jl/stable/ss/#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](https://julianlsolvers.github.io/LsqFit.jl/latest/api/#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. +""" +function system_id_optimal_input_cond( + inp::Union{Vector{Float64}, Matrix{Float64}}, + out::Union{Vector{Float64}, Matrix{Float64}}, + tt::Vector{Float64}, + order::Int, + input_cond::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} + key_list = collect(keys(input_cond_args_guess)) + function model(inp, param) + input_cond_kwargs = Dict(key => param[ii] for (ii, key) ∈ enumerate(key_list)) + sys = system_id( + inp, + out, + tt, + order; + input_cond, + inp_offset, + inp_factor, + out_offset, + out_factor, + input_cond_kwargs, + newpem_kwargs, + verbose, + ) + return model_evolve( + sys, + inp; + input_cond, + inp_offset, + inp_factor, + out_offset, + out_factor, + ) + end + guess = [input_cond_args_guess[key] for key ∈ key_list] + lower = [ + input_cond_args_lower[key] for + key ∈ key_list if key in keys(input_cond_args_lower) + ] + lower = [ + input_cond_args_upper[key] for + key ∈ key_list if key in keys(input_cond_args_upper) + ] + + fit_result = coef(curve_fit(model, inp, out, guess)) + opt = Dict{Symbol, T}(key => fit_result[ii] for (ii, key) ∈ enumerate(key_list)) + + sys = system_id( + inp, + out, + tt, + order; + input_cond, + inp_offset, + inp_factor, + out_offset, + out_factor, + input_cond_kwargs=opt, + newpem_kwargs, + verbose, + ) + return sys, opt +end diff --git a/test/runtests.jl b/test/runtests.jl index 65bfdb0..6060d7d 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,3 +1,4 @@ +using SOLPS2ctrl using SOLPS2ctrl: SOLPS2ctrl using SOLPS2imas: SOLPS2imas using IMAS: IMAS @@ -8,6 +9,8 @@ using Unitful: Unitful using Interpolations: Interpolations using ArgParse: ArgParse using IMASggd: IMASggd, get_grid_subset +using DelimitedFiles: readdlm +using Statistics: mean function parse_commandline() # Define newARGS = ["--yourflag"] to run only tests on your flags when including runtests.jl @@ -36,6 +39,9 @@ function parse_commandline() ["--preparation"], Dict(:help => "Test only preparation", :action => :store_true), + ["--system_id"], + Dict(:help => "Test only system identification", + :action => :store_true), ) args = ArgParse.parse_args(localARGS, s) if !any(values(args)) # If no flags are set, run all tests @@ -305,7 +311,10 @@ if args["heavy_utilities"] dd.core_profiles.profiles_1d[1].electrons.density, dd.core_profiles.profiles_1d[1], dd.equilibrium.time_slice[1], - ).(r, z) + ).( + r, + z, + ) @test size(density_on_grid) == (length(rg), length(zg)) end end @@ -459,3 +468,135 @@ if args["preparation"] @time IMAS.imas2json(dd, filename * ".json", strict=true, freeze=false) end end + +if args["system_id"] + @testset "system_id" begin + case = "$(@__DIR__)/../sample/D3D_Lore_Time_Dep" + println("Reading interferometer data...") + print("json2imas time: ") + @time ids = IMAS.json2imas(case * "/imas_with_ifo.json") + # The file above "imas_with_ifo.json" has been created using following code. + # We jump directly to after the interferometer addition to avoid making + # FusionSyntheticDiagnostics a dependency of this repo. + # + # data_dir = case * "/puff2.5e21_ss_noLat_dribble_308_2d_output" + # b2fgmtry = case * "/baserun/b2fgmtry" + # b2time = data_dir * "/b2time.nc" + # b2mn = data_dir * "/b2mn.dat" + # fort = (data_dir * "/fort.33", data_dir * "/fort.34", data_dir * "/fort.35") + # eqdsk = case * "/baserun/gfile" + # SOLPS2ctrl.geqdsk_to_imas!(eqdsk, ids) + # ids = SOLPS2imas.solps2imas(b2fgmtry, b2time; b2mn=b2mn, fort=fort) + # SOLPS2ctrl.fill_in_extrapolated_core_profile!( + # ids, + # "electrons.density"; + # method="simple", + # cell_subset_idx=-5, + # ) + # FusionSyntheticDiagnostics.add_interferometer!( + # case * "/interferometer.json", + # ids; + # n_e_gsi=-5, + # overwrite=true, + # ) + + ifo_ch = 3 + ne_factor = 1e-19 + ne_offset = ids.interferometer.channel[ifo_ch].n_e_line_average.data[1] + boltzmann_k = 1.38064852e-23 + Temp = 300.0 + + ifo_tt = ids.interferometer.channel[1].n_e_line_average.time + # Make a uniform linear time series + tt = collect(ifo_tt[1]:0.001:ifo_tt[end]) + + # Input gas data + input_gas_data = readdlm(case * "/input_gas_atomsps.txt"; comments=true) + gas_factor = 1e-21 + gas_offset = input_gas_data[1, 2] + # Intepolate input data along the uniform time + input_gas = + Interpolations.linear_interpolation( + input_gas_data[:, 1] .- 5.0, + input_gas_data[:, 2], + ).( + tt, + ) + + # Output interferometer data + ne_ifo = ids.interferometer.channel[ifo_ch].n_e_line_average.data + # Interpolate output data along the uniform time + ne_uniform = Interpolations.linear_interpolation(ifo_tt, ne_ifo).(tt) + + function inp_cond_step( + gas_so::Float64, + previous_gas_so::Float64, + corr::Float64; + p::Union{Float64, Vector{Float64}}=-0.17, + )::Tuple{Float64, Float64} + if p isa Array + param = p[1] + else + param = p + end + if gas_so + corr > previous_gas_so + corr = corr + (gas_so + corr - previous_gas_so) * param + end + return gas_so + corr, corr + end + + """ + inp_cond(inp_gas; p=0.242) + + Condition the input to have this non-linearity. Only when the gas is increasing, the + increase amount is changed by a factor of 'p' of the actual change. This gives a + response of gas input when it is rising vs falling and thus constitute a non-linearity. + """ + function inp_cond( + gas_so::Array{Float64}; + p::Union{Float64, Vector{Float64}}=-0.17, + )::Array{Float64} + ret_vec = deepcopy(gas_so) + corr = 0.0 + for ii ∈ 2:length(ret_vec) + ret_vec[ii], corr = inp_cond_step(ret_vec[ii], ret_vec[ii-1], corr; p) + end + return ret_vec + end + + # verbose only during CI testinf + verbose = get(ENV, "GITHUB_ACTIONS", "") == "true" + + # Third order linear fit of the model + linear_plant_3 = SOLPS2ctrl.system_id( + input_gas, ne_uniform, tt, 3; + inp_offset=gas_offset, inp_factor=gas_factor, + out_offset=ne_offset, out_factor=ne_factor, + verbose, + ) + + # Non-linear input + 3rd order linear fit + nonlin_plant_3, p_opt = system_id_optimal_input_cond( + input_gas, ne_uniform, tt, 3, inp_cond, Dict{Symbol, Any}(:p => -0.2); + inp_offset=gas_offset, inp_factor=gas_factor, + out_offset=ne_offset, out_factor=ne_factor, + verbose, + ) + + lin_out = model_evolve( + linear_plant_3, input_gas; + inp_offset=gas_offset, inp_factor=gas_factor, + out_offset=ne_offset, out_factor=ne_factor, + ) + + nonlin_out = model_evolve( + nonlin_plant_3, input_gas; + inp_offset=gas_offset, inp_factor=gas_factor, + out_offset=ne_offset, out_factor=ne_factor, + input_cond=inp_cond, input_cond_kwargs=p_opt) + + @test sqrt(mean((lin_out - ne_uniform) .^ 2)) < 1.2e18 + + @test sqrt(mean((nonlin_out - ne_uniform) .^ 2)) < 0.4e18 + end +end