Adding System Identification Tools

Project.toml

@@ -1,16 +1,21 @@
 name = "SOLPS2ctrl"
 uuid = "a531d12f-ac8a-43e8-b6d9-bd121431dd49"
 authors = ["David Eldon <eldond@fusion.gat.com>"]
-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"

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

sample/D3D_Lore_Time_Dep/.gitignore

@@ -0,0 +1,2 @@
+/imas_with_ifo.json
+/input_gas_atomsps.txt

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

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

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)

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)

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] = q_extend[rho_output.<rf]
-    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