Replacing OMAS with IMASDD everywhere

Project.toml

@@ -8,6 +8,7 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63"
 Contour = "d38c429a-6771-53c6-b99e-75d170b6e991"
 EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17"
 GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def"
+IMASDD = "06b86afa-9f21-11ec-2ef8-e51b8960cfc5"
 Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f"

makefile

@@ -15,16 +15,17 @@ env_with_cloned_repo r:
 	-dn=$(shell dirname $(shell pwd)); \
 	if [[ "$${dn:(-10)}" == ".julia/dev" ]]; then ext="" ; else ext=".jl";fi; \
 	git clone "git@github.com:ProjectTorreyPines/OMAS.jl.git" ../OMAS$${ext}; \
+	git clone "git@github.com:ProjectTorreyPines/IMASDD.jl.git" ../IMASDD$${ext}; \
 	git clone "git@github.com:ProjectTorreyPines/GGDUtils.jl.git" ../GGDUtils$${ext}; \
 	git clone "git@github.com:ProjectTorreyPines/SOLPS2IMAS.jl.git" ../SOLPS2IMAS$${ext}; \
-	julia --project=. -e 'using Pkg; Pkg.rm(["OMAS", "GGDUtils", "SOLPS2IMAS", "EFIT"]); Pkg.develop(path="../OMAS'$${ext}'"); Pkg.develop(path="../GGDUtils'$${ext}'"); Pkg.develop(path="../SOLPS2IMAS'$${ext}'"); Pkg.add(url="git@github.com:JuliaFusion/EFIT.jl.git", rev="master"); Pkg.instantiate()'
+	julia --project=. -e 'using Pkg; Pkg.rm(["OMAS", "IMASDD", "GGDUtils", "SOLPS2IMAS", "EFIT"]); Pkg.develop(path="../OMAS'$${ext}'"); Pkg.develop(path="../IMASDD'$${ext}'"); Pkg.develop(path="../GGDUtils'$${ext}'"); Pkg.develop(path="../SOLPS2IMAS'$${ext}'"); Pkg.add(url="git@github.com:JuliaFusion/EFIT.jl.git", rev="master"); Pkg.instantiate()'
 
 env_with_git_url u:
 	@echo "Pulling sample files using dvc"
 	-dvc pull
 	@echo "Creating Julia environment with the git urls without creating local clones"
 	@echo "Generating Project.toml and Manifest.toml"
-	julia --project=. -e 'using Pkg; Pkg.rm(["OMAS", "GGDUtils", "SOLPS2IMAS", "EFIT"]); Pkg.add(url="git@github.com:ProjectTorreyPines/OMAS.jl.git", rev="master"); Pkg.add(url="git@github.com:ProjectTorreyPines/GGDUtils.jl.git", rev="master"); Pkg.add(url="git@github.com:ProjectTorreyPines/SOLPS2IMAS.jl.git", rev="master"); Pkg.add(url="git@github.com:JuliaFusion/EFIT.jl.git", rev="master"); Pkg.instantiate()'
+	julia --project=. -e 'using Pkg; Pkg.rm(["OMAS", "IMASDD", "GGDUtils", "SOLPS2IMAS", "EFIT"]); Pkg.add(url="git@github.com:ProjectTorreyPines/OMAS.jl.git", rev="master"); Pkg.add(url="git@github.com:ProjectTorreyPines/IMASDD.jl.git", rev="master"); Pkg.add(url="git@github.com:ProjectTorreyPines/GGDUtils.jl.git", rev="master"); Pkg.add(url="git@github.com:ProjectTorreyPines/SOLPS2IMAS.jl.git", rev="master"); Pkg.add(url="git@github.com:JuliaFusion/EFIT.jl.git", rev="master"); Pkg.instantiate()'
 
 clean:
 	@echo "Deleting Manifest.toml"

sample/geqdsk_iter_small_sample

@@ -1,4 +1,4 @@
-  EFITD    00/00/2008    #002296  0200ms Convert   0  65  65
+  EFITD    00/00/2008    #002296  0200   0  65  65
  0.624390220E+01 0.124878049E+02 0.619999981E+01 0.307804870E+01 0.200000050E+00
  0.635000000E+01 0.580000000E+00 0.118000000E+02 0.000000000E+00 0.530000019E+01
  0.150229038E+08 0.118000000E+02 0.000000000E+00 0.635000000E+01 0.000000000E+00

src/SD4SOLPS.jl

@@ -1,6 +1,6 @@
 module SD4SOLPS
 
-using OMAS: OMAS
+import OMAS as IMASDD
 using SOLPS2IMAS: SOLPS2IMAS
 using EFIT: EFIT
 using Interpolations: Interpolations
@@ -68,13 +68,20 @@ end
 Transfers the equilibrium reconstruction in an EFIT-style gEQDSK file into
 the IMAS DD structure.
 """
-function geqdsk_to_imas!(eqdsk_file, dd; time_index=1)
+function geqdsk_to_imas!(eqdsk_file, dd; set_time=nothing, time_index=1)
     # https://github.com/JuliaFusion/EFIT.jl/blob/master/src/io.jl
-    g = EFIT.readg(eqdsk_file)
+    g = EFIT.readg(eqdsk_file; set_time=set_time)
     # Copying ideas from OMFIT: omfit/omfit_classes/omfit_eqdsk.py / to_omas()
     eq = dd.equilibrium
-    resize!(eq.time_slice, 1)
+    if IMASDD.ismissing(eq, :time)
+        eq.time = Array{Float64}(undef, time_index)
+    end
+    eq.time[time_index] = g.time
+    if length(eq.time_slice) < time_index
+        resize!(eq.time_slice, time_index)
+    end
     eqt = eq.time_slice[time_index]
+    eqt.time = g.time
 
     # 0D
     gq = eqt.global_quantities
@@ -287,7 +294,7 @@ function preparation(
     println("Extrapolated edge profiles (but not really (placeholder only))")
 
     if output_format == "json"
-        OMAS.imas2json(dd, filename * ".json")
+        IMASDD.imas2json(dd, filename * ".json")
     else
         throw(ArgumentError(string("Unrecognized output format: ", output_format)))
     end

src/repair_eq.jl

@@ -17,7 +17,7 @@ export check_rho_1d
 
 Checks to see if rho exists and is valid in the equilibrium 1d profiles
 """
-function check_rho_1d(dd::OMAS.dd; time_slice::Int64=1, throw_on_fail::Bool=false)
+function check_rho_1d(dd::IMASDD.dd; time_slice::Int64=1, throw_on_fail::Bool=false)
     rho = dd.equilibrium.time_slice[time_slice].profiles_1d.rho_tor_norm
     if length(rho) < 1
         rho_okay = false
@@ -48,11 +48,11 @@ function check_rho_1d(dd::OMAS.dd; time_slice::Int64=1, throw_on_fail::Bool=fals
 end
 
 """
-    function add_rho_to_equilibrium(dd:OMAS.dd)
+    function add_rho_to_equilibrium(dd:IMASDD.dd)
 
 Adds equilibrium rho profile to the DD
 """
-function add_rho_to_equilibrium!(dd::OMAS.dd)
+function add_rho_to_equilibrium!(dd::IMASDD.dd)
     nt = length(dd.equilibrium.time_slice)
     if nt < 1
         println("No equilibrium time slices to work with; can't add rho")
@@ -71,8 +71,14 @@ function add_rho_to_equilibrium!(dd::OMAS.dd)
                 continue
             end
         end
-        if length(eqt.profiles_1d.phi) == 0
-            resize!(eqt.profiles_1d.phi, n)
+        if (
+            if IMASDD.ismissing(eqt.profiles_1d, :phi)
+                true
+            else
+                IMASDD.isempty(eqt.profiles_1d.phi)
+            end
+        )
+            eqt.profiles_1d.phi = Array{Float64}(undef, n)
             psi2 = eqt.profiles_2d[1].psi
             req = collect(eqt.profiles_2d[1].grid.dim1)
             zeq = collect(eqt.profiles_2d[1].grid.dim2)

src/supersize_profile.jl

@@ -3,7 +3,7 @@ Utilities for extrapolating profiles
 """
 
 # import CalculusWithJulia
-using OMAS: OMAS
+import OMAS as IMASDD
 using Interpolations: Interpolations
 using GGDUtils:
     GGDUtils, get_grid_subset, add_subset_element!, get_subset_boundary,
@@ -55,7 +55,7 @@ Concept:
     integrate it to get the profile of the quantity in question
 """
 function extrapolate_core(edge_rho, edge_quantity, rho_output)
-    grad = OMAS.gradient(edge_rho, edge_quantity)
+    grad = IMASDD.gradient(edge_rho, edge_quantity)
     gf = grad[1]
     rf = edge_rho[1]
     gmid = -abs(gf) / 4.0
@@ -92,7 +92,7 @@ end
 #!format off
 """
     fill_in_extrapolated_core_profile!(
-    dd::OMAS.dd,
+    dd::IMASDD.dd,
     quantity_name::String;
     method::String="simple",
     eq_time_idx::Int64=1,
@@ -108,7 +108,7 @@ edge_profiles data to rho, calls the function that performs the extrapolation (w
 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.
-dd: an IMAS/OMAS data dictionary
+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.
@@ -130,7 +130,7 @@ cell_subset_idx: index of the subset of cells to use for the extrapolation. The
 """
 #!format on
 function fill_in_extrapolated_core_profile!(
-    dd::OMAS.dd,
+    dd::IMASDD.dd,
     quantity_name::String;
     method::String="simple",
     eq_time_idx::Int64=1,
@@ -231,10 +231,9 @@ function fill_in_extrapolated_core_profile!(
             ).(psi_for_quantity[in_bounds])
 
         # Make sure the output 1D rho grid exists; create it if needed
-        if length(dd.core_profiles.profiles_1d[it].grid.rho_tor_norm) == 0
-            resize!(dd.core_profiles.profiles_1d[it].grid.rho_tor_norm, 201)
-            # If you don't like this default, then you should write grid.rho_tor_norm before
-            # calling this function.
+        if IMASDD.ismissing(dd.core_profiles.profiles_1d[it].grid, :rho_tor_norm)
+            # If you don't like this default, then you should write grid.rho_tor_norm
+            # before calling this function.
             dd.core_profiles.profiles_1d[it].grid.rho_tor_norm =
                 collect(LinRange(0, 1, 201))
         end
@@ -300,7 +299,7 @@ Returns:
   - normalized poloidal flux on the equilibrium grid
   - a linear interpolation of norm pol flux vs. R and Z, ready to be evaluated
 """
-function prep_flux_map(dd::OMAS.dd; eq_time_idx::Int64=1, eq_profiles_2d_idx::Int64=1)
+function prep_flux_map(dd::IMASDD.dd; eq_time_idx::Int64=1, eq_profiles_2d_idx::Int64=1)
     eqt = dd.equilibrium.time_slice[eq_time_idx]
     p2 = eqt.profiles_2d[eq_profiles_2d_idx]
     r_eq = p2.grid.dim1
@@ -316,7 +315,7 @@ end
 #! format off
 """
     mesh_psi_spacing(
-    dd::OMAS.dd;
+    dd::IMASDD.dd;
     eq_time_idx::Int64=1,
     eq_profiles_2d_idx::Int64=1,
     grid_ggd_idx::Int64=1,
@@ -346,7 +345,7 @@ spacing_rule: "edge" or "mean" to make spacing of new cells (in psi_N) be the sa
 """
 #! format on
 function mesh_psi_spacing(
-    dd::OMAS.dd;
+    dd::IMASDD.dd;
     eq_time_idx::Int64=1,
     eq_profiles_2d_idx::Int64=1,
     grid_ggd_idx::Int64=1,
@@ -417,7 +416,7 @@ out to the most distant (in flux space) point on the limiting surface.
 Returns a vector of psi_N levels.
 """
 function pick_extension_psi_range(
-    dd::OMAS.dd;
+    dd::IMASDD.dd;
     eq_time_idx::Int64=1,
     eq_profiles_2d_idx::Int64=1,
     grid_ggd_idx::Int64=1,
@@ -553,7 +552,7 @@ function mesh_ext_follow_grad(
     end
 
     # Step along the paths of steepest descent to populate the mesh.
-    dpsindr, dpsindz = OMAS.gradient(r_eq, z_eq, psin_eq)
+    dpsindr, dpsindz = IMASDD.gradient(r_eq, z_eq, psin_eq)
     dpdr = Interpolations.linear_interpolation((r_eq, z_eq), dpsindr)
     dpdz = Interpolations.linear_interpolation((r_eq, z_eq), dpsindz)
     rlim = (minimum(r_eq), maximum(r_eq))
@@ -596,7 +595,7 @@ mesh_r: matrix of R values for the extended mesh
 mesh_z: matrix of Z values for the extended mesh
 """
 function modify_mesh_ext_near_x!(
-    eqt::OMAS.equilibrium__time_slice,
+    eqt::IMASDD.equilibrium__time_slice,
     mesh_r::Matrix{Float64},
     mesh_z::Matrix{Float64},
 )
@@ -750,15 +749,9 @@ function record_regular_mesh!(
 
     # Preserve record of standard (non extended) mesh
     for i ∈ 1:5
-        std_sub = get_grid_subset(grid_ggd, -i)
         orig_sub = get_grid_subset(grid_ggd, i)
-        resize!(std_sub.element, length(orig_sub.element))
-        for j ∈ 1:length(orig_sub.element)
-            std_sub.element[j] = deepcopy(orig_sub.element[j])
-        end
-        std_sub.identifier.index = -i
-        std_sub.dimension = deepcopy(orig_sub.dimension)
-        std_sub.metric = deepcopy(orig_sub.metric)
+        grid_ggd.grid_subset[n_existing_subsets+i] = deepcopy(orig_sub)
+        grid_ggd.grid_subset[n_existing_subsets+i].identifier.index = -i
     end
     all_nodes_sub = get_grid_subset(grid_ggd, 1)
     all_edges_sub = get_grid_subset(grid_ggd, 2)
@@ -873,7 +866,7 @@ clear_cache: delete any existing cache file (for use in testing)
 """
 #!format on
 function cached_mesh_extension!(
-    dd::OMAS.dd,
+    dd::IMASDD.dd,
     eqdsk_file::String,
     b2fgmtry::String;
     eq_time_idx::Int64=1,
@@ -942,7 +935,7 @@ end
 Extends the mesh out into the SOL
 """
 function mesh_extension_sol!(
-    dd::OMAS.dd;
+    dd::IMASDD.dd;
     eq_time_idx::Int64=1,
     eq_profiles_2d_idx::Int64=1,
     grid_ggd_idx::Int64=1,
@@ -989,13 +982,13 @@ end
 
 """
     fill_in_extrapolated_edge_profile!(
-        dd::OMAS.dd, quantity_name::String; method::String="simple",
+        dd::IMASDD.dd, quantity_name::String; method::String="simple",
     )
 
 JUST A PLACEHOLDER FOR NOW. DOESN'T ACTUALLY WORK YET.
 """
 function fill_in_extrapolated_edge_profile!(
-    dd::OMAS.dd,
+    dd::IMASDD.dd,
     quantity_name::String;
     method::String="simple",
     eq_time_idx::Int64=1,

test/runtests.jl

@@ -1,6 +1,6 @@
 using SD4SOLPS: SD4SOLPS
 using SOLPS2IMAS: SOLPS2IMAS
-using OMAS: OMAS
+import OMAS as IMASDD
 using EFIT: EFIT
 using Plots
 using Test
@@ -183,8 +183,9 @@ if args["core_profile_extension"]
         @test isfile(b2time)
         @test isfile(b2mn)
         @test isfile(eqdsk)
+        eqdsk_time = parse(Float64, split(eqdsk, ".")[end]) / 1000.0
         dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time; b2mn=b2mn)
-        SD4SOLPS.geqdsk_to_imas!(eqdsk, dd)
+        SD4SOLPS.geqdsk_to_imas!(eqdsk, dd; set_time=eqdsk_time)
         rho = dd.equilibrium.time_slice[1].profiles_1d.rho_tor_norm
 
         if !SD4SOLPS.check_rho_1d(dd; time_slice=1)
@@ -219,8 +220,9 @@ if args["edge_profile_extension"]
     @testset "edge_profile_extension" begin
         # Test for getting mesh spacing
         b2fgmtry, b2time, b2mn, eqdsk = define_default_sample_set()
+        eqdsk_time = parse(Float64, split(eqdsk, ".")[end]) / 1000.0
         dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time; b2mn=b2mn)
-        SD4SOLPS.geqdsk_to_imas!(eqdsk, dd)
+        SD4SOLPS.geqdsk_to_imas!(eqdsk, dd; set_time=eqdsk_time)
         dpsin = SD4SOLPS.mesh_psi_spacing(dd)
         @test dpsin > 0.0
 
@@ -289,7 +291,7 @@ if args["heavy_utilities"]
 
         # Test for sweeping 1D core profiles into 2D R,Z
         # (or anyway evaluating them at any R,Z location)
-        dd = OMAS.dd()
+        dd = IMASDD.dd()
         eqdsk_file =
             splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample"
         SD4SOLPS.geqdsk_to_imas!(eqdsk_file, dd)
@@ -298,8 +300,6 @@ if args["heavy_utilities"]
         resize!(dd.core_profiles.profiles_1d, prof_time_idx)
         n = 101
         rho_n = Array(LinRange(0, 1.0, n))
-        resize!(dd.core_profiles.profiles_1d[prof_time_idx].grid.rho_tor_norm, n)
-        resize!(dd.core_profiles.profiles_1d[prof_time_idx].electrons.density, n)
         dd.core_profiles.profiles_1d[prof_time_idx].grid.rho_tor_norm = rho_n
         dd.core_profiles.profiles_1d[prof_time_idx].electrons.density =
             make_test_profile(rho_n)
@@ -323,7 +323,7 @@ end
 if args["repair_eq"]
     @testset "repair_eq" begin
         # Prepare sample
-        dd = OMAS.dd()
+        dd = IMASDD.dd()
         eqdsk = splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample"
         SD4SOLPS.geqdsk_to_imas!(eqdsk, dd)
         # Make sure rho is missing
@@ -355,8 +355,18 @@ if args["geqdsk_to_imas"]
         tslice = 1
         for sample_file ∈ sample_files
             println(sample_file)
-            dd = OMAS.dd()
-            SD4SOLPS.geqdsk_to_imas!(sample_file, dd; time_index=tslice)
+            dd = IMASDD.dd()
+            if endswith(sample_file, "00")
+                eqdsk_time = parse(Float64, split(sample_file, ".")[end]) / 1000.0
+            else
+                eqdsk_time = nothing
+            end
+            SD4SOLPS.geqdsk_to_imas!(
+                sample_file,
+                dd;
+                set_time=eqdsk_time,
+                time_index=tslice,
+            )
             eqt = dd.equilibrium.time_slice[tslice]
 
             # global