diff --git a/Project.toml b/Project.toml index a3e1890..7c1d4a6 100644 --- a/Project.toml +++ b/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" diff --git a/makefile b/makefile index 53b0b8f..197d93b 100644 --- a/makefile +++ b/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" diff --git a/sample/geqdsk_iter_small_sample b/sample/geqdsk_iter_small_sample index e51be73..a230a97 100644 --- a/sample/geqdsk_iter_small_sample +++ b/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 diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index ae713b3..75bfe79 100644 --- a/src/SD4SOLPS.jl +++ b/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 diff --git a/src/repair_eq.jl b/src/repair_eq.jl index 54e7937..48e2100 100644 --- a/src/repair_eq.jl +++ b/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) diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index b541c75..e1c6a68 100644 --- a/src/supersize_profile.jl +++ b/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, diff --git a/test/runtests.jl b/test/runtests.jl index e2551c1..70fd0ed 100644 --- a/test/runtests.jl +++ b/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