diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml new file mode 100644 index 0000000..aa7f0eb --- /dev/null +++ b/.JuliaFormatter.toml @@ -0,0 +1,27 @@ +align_assignment = true +align_conditional = true +align_matrix = true +align_pair_arrow = true +align_struct_field = true +always_for_in = true +always_use_return = true +annotate_untyped_fields_with_any = false +conditional_to_if = true +for_in_replacement = "∈" +format_docstrings = true +import_to_using = true +indent = 4 +indent_submodule = true +join_lines_based_on_source = true +long_to_short_function_def = true +margin = 88 +normalize_line_endings = "unix" +pipe_to_function_call = true +remove_extra_newlines = true +separate_kwargs_with_semicolon = true +surround_whereop_typeparameters = true +trailing_comma = true +whitespace_in_kwargs = false +whitespace_ops_in_indices = false +whitespace_typedefs = true +yas_style_nesting = true diff --git a/.github/workflows/format_check.yml b/.github/workflows/format_check.yml new file mode 100644 index 0000000..4e46a44 --- /dev/null +++ b/.github/workflows/format_check.yml @@ -0,0 +1,36 @@ +name: Format Check + +on: + push: + branches: ["master", "dev", "format"] + pull_request: + branches: ["master", "dev"] +jobs: + check: + runs-on: ${{ matrix.os }} + strategy: + matrix: + julia-version: [1.9.3] + julia-arch: [x86] + os: [ubuntu-latest] + steps: + - uses: julia-actions/setup-julia@latest + with: + version: ${{ matrix.julia-version }} + + - uses: actions/checkout@v1 + - name: Install JuliaFormatter and format + run: | + julia -e 'using Pkg; Pkg.add(PackageSpec(name="JuliaFormatter"))' + julia -e 'using JuliaFormatter; format(".", verbose=true)' + - name: Format check + run: | + julia -e ' + out = Cmd(`git diff --name-only`) |> read |> String + if out == "" + exit(0) + else + @error "Some files have not been formatted !!!" + write(stdout, out) + exit(1) + end' diff --git a/.gitignore b/.gitignore index 6b72f03..73c794c 100644 --- a/.gitignore +++ b/.gitignore @@ -3,3 +3,4 @@ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml +sd_input_data.json diff --git a/sample/ITER_Lore_2296_00000.dvc b/sample/ITER_Lore_2296_00000.dvc index 01475a8..6c44f0d 100644 --- a/sample/ITER_Lore_2296_00000.dvc +++ b/sample/ITER_Lore_2296_00000.dvc @@ -1,10 +1,10 @@ -md5: fd0856e1e06a7ddab0b8e4aaddc170de +md5: 377b6f3a1439c78e6a310739b4ca69d0 frozen: true deps: - path: ITER_Lore_2296_00000 repo: - url: git@github.com:ProjectTorreyPines/TestSamples.git - rev_lock: f905bd19e48bd54bc53416de9bbdaadda6654b9a + url: git@github.com:ProjectTorreyPines/SOLPSTestSamples.git + rev_lock: 35fd24298a4118eec6491dcd17d4e0ec97823e83 outs: - md5: 48af1b920caaa6bc1788af8a36512638.dir size: 322802293 diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 089c493..6f46882 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -1,9 +1,9 @@ module SD4SOLPS -import OMAS -import SOLPS2IMAS -import EFIT -import Interpolations +using OMAS: OMAS +using SOLPS2IMAS: SOLPS2IMAS +using EFIT: EFIT +using Interpolations: Interpolations #import GGDUtils export find_files_in_allowed_folders @@ -21,22 +21,32 @@ Searches a list of allowed folders for a set of filenames that will provide info about the SOLPS case. Returns a list of filenames with complete paths. Example: -SD4SOLPS.find_files_in_allowed_folders("/D3D_Ma_184833_03600", eqdsk_file="g184833.03600") +SD4SOLPS.find_files_in_allowed_folders( +"/D3D_Ma_184833_03600", +eqdsk_file="g184833.03600", +) """ -function find_files_in_allowed_folders(input_dirs...; eqdsk_file, recursive=true, allow_reduced_versions=false) +function find_files_in_allowed_folders( + input_dirs...; + eqdsk_file, + recursive=true, + allow_reduced_versions=false, +) files = ["b2fgmtry", "b2time.nc", "b2mn.dat", "gridspacedesc.yml", eqdsk_file] - reduced_files = ["b2fgmtry_red", "b2time_red.nc", "b2mn.dat", "gridspacedesc.yml", eqdsk_file] + reduced_files = + ["b2fgmtry_red", "b2time_red.nc", "b2mn.dat", "gridspacedesc.yml", eqdsk_file] output_files = fill("", length(files)) if recursive dirs = [] - for dir in input_dirs - dirs = append!(dirs, [subdir[1] for subdir in [item for item in walkdir(dir)]]) + for dir ∈ input_dirs + dirs = + append!(dirs, [subdir[1] for subdir ∈ [item for item ∈ walkdir(dir)]]) end else dirs = input_dirs end - for i in 1:length(files) - for dir in dirs + for i ∈ 1:length(files) + for dir ∈ dirs file = dir * "/" * files[i] reduced_file = dir * "/" * reduced_files[i] if isfile(file) @@ -87,7 +97,8 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1) p1.f_df_dpsi = g.ffprim p1.dpressure_dpsi = g.pprime p1.q = g.qpsi - if hasproperty(g, :rhovn) # rhovn is not in the original EFIT.jl but is added on a branch + if hasproperty(g, :rhovn) + # rhovn is not in the original EFIT.jl but is added on a branch p1.rho_tor_norm = g.rhovn end @@ -122,8 +133,7 @@ function geqdsk_to_imas(eqdsk_file, dd; time_index=1) limiter.type.description = "first wall" resize!(limiter.unit, 1) limiter.unit[1].outline.r = g.rlim - limiter.unit[1].outline.z = g.zlim - + return limiter.unit[1].outline.z = g.zlim end """ @@ -136,14 +146,14 @@ r: R locations of desired output points / m z: Z locations of desired output points / m """ function core_profile_2d(dd, prof_time_idx, eq_time_idx, quantity, r, z) - if !check_rho_1d(dd, time_slice=eq_time_idx) + if !check_rho_1d(dd; time_slice=eq_time_idx) throw(ArgumentError("Equilibrium rho profile in input DD was missing.")) end prof = dd.core_profiles.profiles_1d[prof_time_idx] rho_prof = prof.grid.rho_tor_norm quantity_fields = split(quantity, ".") p = prof - for field in quantity_fields + for field ∈ quantity_fields p = getproperty(p, Symbol(field)) end eqt = dd.equilibrium.time_slice[eq_time_idx] @@ -182,14 +192,16 @@ function core_profile_2d(dd, prof_time_idx, eq_time_idx, quantity, r, z) # EQUILBRIUM (the connection from rho to R,Z via psi): # psin_eq_ext : 1D array of psi_N values that corresponds to rhon_eq_ext - # rhon_eq_ext : 1D array of rho_N values that corresponds to psin_eq_ext --> connects rho to everything else via psi + # rhon_eq_ext : 1D array of rho_N values that corresponds to psin_eq_ext + # --> connects rho to everything else via psi # psinrz : 2D array of psi_N values that corresponds to r_eq and z_eq # r_eq and z_eq : 1D coordinate arrays that correspond to psinrz # OUTPUT INSTRUCTIONS: # r and z : coordinates of output points where values of p are desired - psi_at_requested_points = Interpolations.LinearInterpolation((r_eq, z_eq), psinrz).(r, z) + psi_at_requested_points = + Interpolations.LinearInterpolation((r_eq, z_eq), psinrz).(r, z) rhonpsi = Interpolations.LinearInterpolation(psin_eq_ext, rhon_eq_ext) rho_at_requested_points = rhonpsi.(psi_at_requested_points) itp = Interpolations.LinearInterpolation(rho_prof, p) @@ -203,19 +215,41 @@ end Gathers SOLPS and EFIT files and loads them into IMAS structure. Extrapolates profiles as needed to get a complete picture. """ -function preparation(eqdsk_file, dirs...) - b2fgmtry, b2time, b2mn, gridspec, eqdsk = find_files_in_allowed_folders( - dirs..., eqdsk_file=eqdsk_file - ) +function preparation( + eqdsk_file, + dirs...; + core_method::String="simple", + edge_method::String="simple", + filename::String="sd_input_data", + output_format::String="json", +) + b2fgmtry, b2time, b2mn, gridspec, eqdsk = + find_files_in_allowed_folders(dirs...; eqdsk_file=eqdsk_file) println("Found source files:") println(" b2fgmtry = ", b2fgmtry) println(" b2time = ", b2time) println(" b2mn.dat = ", b2mn) println(" gridspec = ", gridspec) println(" eqdsk = ", eqdsk) + dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) geqdsk_to_imas(eqdsk, dd) - println("Loaded data into IMAS DD") + println("Loaded input data into IMAS DD") + + fill_in_extrapolated_core_profile!(dd, "electrons.density"; method=core_method) + fill_in_extrapolated_core_profile!(dd, "electrons.temperature"; method=core_method) + # ... more profiles here as they become available in b2time + println("Extrapolated core profiles") + + fill_in_extrapolated_edge_profile!(dd, "electrons.density"; method=core_method) + # ... more profiles here + println("Extrapolated edge profiles (but not really (placeholder only))") + + if output_format == "json" + OMAS.imas2json(dd, filename * ".json") + else + throw(ArgumentError(string("Unrecognized output format: ", output_format))) + end return dd end diff --git a/src/repair_eq.jl b/src/repair_eq.jl index 08feddb..7954c4f 100644 --- a/src/repair_eq.jl +++ b/src/repair_eq.jl @@ -6,8 +6,8 @@ use this tool mainly for test cases, as a utility. For a real equilibrium, problems should be fixed properly. """ -import Contour -import Statistics +using Contour: Contour +using Statistics: Statistics export add_rho_to_equilibrium! export check_rho_1d @@ -40,7 +40,7 @@ function add_rho_to_equilibrium!(dd::OMAS.dd) println("No equilibrium time slices to work with; can't add rho") return end - for it = 1:nt + for it ∈ 1:nt eqt = dd.equilibrium.time_slice[it] b0 = dd.equilibrium.vacuum_toroidal_field.b0[it] r0 = dd.equilibrium.vacuum_toroidal_field.r0[it] @@ -60,7 +60,9 @@ function add_rho_to_equilibrium!(dd::OMAS.dd) zeq = collect(eqt.profiles_2d[1].grid.dim2) if (length(req), length(zeq)) == size(psi2') psi2 = Matrix(psi2') - println("transposed psi to make it compatible with r,z prior to contouring") + println( + "transposed psi to make it compatible with r,z prior to contouring", + ) end if (length(req), length(zeq)) != size(psi2) println("Invalid equilibrium data. rho cannot be added.") @@ -70,25 +72,38 @@ function add_rho_to_equilibrium!(dd::OMAS.dd) println(" psi2d : ", size(psi2)) return else - println("Eq looks okay ", (length(req), length(zeq)), ", ", size(psi2), ". ", (size(req), size(zeq))) + println( + "Eq looks okay ", + (length(req), length(zeq)), + ", ", + size(psi2), + ". ", + (size(req), size(zeq)), + ) end - for j = 1:n + for j ∈ 1:n contour_level = psi[j] if j == n - # The last contour has a X point if everything is lined up right, and that could get - # weird. We want a contour that only takes the core boundary part of the separatrix. + # The last contour has a X point if everything is lined up right, + # and that could get weird. We want a contour that only takes the + # core boundary part of the separatrix. r = eqt.boundary.outline.r z = eqt.boundary.outline.z else c = Contour.contour(req, zeq, psi2, contour_level) clines = Contour.lines(c) - line_avg_height = [Statistics.mean([clines[i].vertices[v][2] for v in 1:length(clines[i].vertices)]) for i in 1:length(clines)] + line_avg_height = [ + Statistics.mean([ + clines[i].vertices[v][2] for + v ∈ 1:length(clines[i].vertices) + ]) for i ∈ 1:length(clines) + ] cl = clines[argmin(abs.(line_avg_height))] # Now just do the 2D integral of B within the area of the contour - r = [cl.vertices[v][1] for v in 1:length(cl.vertices)] - z = [cl.vertices[v][2] for v in 1:length(cl.vertices)] + r = [cl.vertices[v][1] for v ∈ 1:length(cl.vertices)] + z = [cl.vertices[v][2] for v ∈ 1:length(cl.vertices)] end - # Get the upper and lower halves of the path, each going from rmin to rmax. + # Get the upper & lower halves of the path, each going from rmin to rmax rmin = minimum(r) rmax = maximum(r) irmin = argmin(r) @@ -123,13 +138,17 @@ function add_rho_to_equilibrium!(dd::OMAS.dd) z2i = Interpolations.LinearInterpolation(r2, z2) rr = LinRange(rmin, rmax, 101) 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 in 1:length(rc)] + integral_part_ = [ + log(rr[i+1] / rr[i]) * abs(z1i(rc[i]) - z2i(rc[i])) for + i ∈ 1:length(rc) + ] phi_ = r0 * b0 .* integral_part_ eqt.profiles_1d.phi[j] = sum(phi_) end end eqt.profiles_1d.rho_tor = sqrt.(eqt.profiles_1d.phi / (π * b0)) - eqt.profiles_1d.rho_tor_norm = eqt.profiles_1d.rho_tor / eqt.profiles_1d.rho_tor[end] + eqt.profiles_1d.rho_tor_norm = + eqt.profiles_1d.rho_tor / eqt.profiles_1d.rho_tor[end] end - println("Rho has been added to the equilibrium") + return println("Rho has been added to the equilibrium") end diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index 898c339..7b0e100 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -3,10 +3,14 @@ Utilities for extrapolating profiles """ # import CalculusWithJulia -import OMAS -import Interpolations -import NumericalIntegration -import GGDUtils +using OMAS: OMAS +using Interpolations: Interpolations +using NumericalIntegration: NumericalIntegration +using GGDUtils: GGDUtils + +export extrapolate_core +export fill_in_extrapolated_core_profile +export mesh_psi_spacing """ extrapolate_core(edge_rho, edge_quantity, rho_output) @@ -14,14 +18,15 @@ import GGDUtils 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, + + 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 + 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 + 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, + 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 """ function extrapolate_core(edge_rho, edge_quantity, rho_output) @@ -38,16 +43,25 @@ function extrapolate_core(edge_rho, edge_quantity, rho_output) gg = [0, gmid, gped_enforce, gf] rr = [0, rmid, rped_enforce, rf] # https://www.matecdev.com/posts/julia-interpolation.html - itp = Interpolations.LinearInterpolation(rr, gg, extrapolation_bc=Interpolations.Line()) + itp = Interpolations.LinearInterpolation( + rr, + gg; + extrapolation_bc=Interpolations.Line(), + ) #itp = Interpolations.extrapolate(itp, Interpolations.Line()) g = itp(rho_output) q_extend = NumericalIntegration.cumul_integrate(rho_output, g) - q_offset = edge_quantity[1] - Interpolations.LinearInterpolation(rho_output, q_extend)(rf) + q_offset = + edge_quantity[1] - Interpolations.LinearInterpolation(rho_output, q_extend)(rf) 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] = Interpolations.LinearInterpolation(edge_rho, edge_quantity).(rho_output[rho_output .>= rf]) + output_profile[rho_output.=rf] = + Interpolations.LinearInterpolation( + edge_rho, + edge_quantity, + ).(rho_output[rho_output.>=rf]) return output_profile end @@ -55,34 +69,49 @@ end function fill_in_extrapolated_core_profile(dd::OMAS.dd, quantity_name::String) 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. +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. dd: an IMAS/OMAS data dictionary -quantity_name: the name of a quantity in edge_profiles.profiles_2d and core_profiles.profiles_1d, - such as "electrons.density" +quantity_name: the name of a quantity in edge_profiles.profiles_2d and +core_profiles.profiles_1d, such as "electrons.density" """ -function fill_in_extrapolated_core_profile(dd::OMAS.dd, quantity_name::String) - println("hey hey hey") +function fill_in_extrapolated_core_profile!( + dd::OMAS.dd, + quantity_name::String; + method::String="simple", + eq_time_idx::Int64=1, +) ggd_idx = 1 space_number = 1 space = dd.edge_profiles.grid_ggd[ggd_idx].space[space_number] - cell_subset = SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[ggd_idx], 5) - midplane_subset = SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[ggd_idx], 11) - + cell_subset = + SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[ggd_idx], 5) + midplane_subset = + SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[ggd_idx], 11) + if length(midplane_subset.element) < 1 - throw(ArgumentError(string( - "Midplane subset length was ", length(midplane_subset.element), - ". Unacceptable data in data dictionary." - ))) + throw( + ArgumentError( + string( + "Midplane subset length was ", length(midplane_subset.element), + ". Unacceptable data in data dictionary.", + ), + ), + ) end if length(cell_subset.element) < 1 - throw(ArgumentError(string( - "Cell subset length was ", length(cell_subset.element), - ". Unacceptable data in data dictionary." - ))) + throw( + ArgumentError( + string( + "Cell subset length was ", length(cell_subset.element), + ". Unacceptable data in data dictionary.", + ), + ), + ) end nt = 1 @@ -90,46 +119,67 @@ function fill_in_extrapolated_core_profile(dd::OMAS.dd, quantity_name::String) resize!(dd.core_profiles.profiles_1d, nt) end it = 1 - eq_it = 1 # quantity_in_cells = SOLPS2IMAS.val_obj(dd.edge_profiles.ggd[it], quantity_name, ggd_idx) tags = split(quantity_name, ".") quantity_str = dd.edge_profiles.ggd[it] - for tag in tags + for tag ∈ tags quantity_str = getproperty(quantity_str, Symbol(tag)) end # Quantity is in cells now cell_subset_idx = 5 nqv = length(quantity_str[cell_subset_idx].values) if nqv < 1 - println("Quantity ", quantity_name, " has ", nqv, " elements in subset ", cell_subset_idx, ". Unacceptable.") + println( + "Quantity ", + quantity_name, + " has ", + nqv, + " elements in subset ", + cell_subset_idx, + ". Unacceptable.", + ) return end quantity_str[cell_subset_idx].grid_subset_index = cell_subset_idx - midplane_cell_centers, quantity = GGDUtils.project_prop_on_subset!(quantity_str, cell_subset, midplane_subset, space=space) + midplane_cell_centers, quantity = GGDUtils.project_prop_on_subset!( + quantity_str, + cell_subset, + midplane_subset, + space, + ) # Now quantity is at the outboard midplane # Get the rho values to go with the midplane quantity values - # midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) # Not needed; done by project_prop_on_subset - r = [midplane_cell_centers[i][1] for i in 1:length(midplane_cell_centers)] - z = [midplane_cell_centers[i][2] for i in 1:length(midplane_cell_centers)] + r = [midplane_cell_centers[i][1] for i ∈ 1:length(midplane_cell_centers)] + z = [midplane_cell_centers[i][2] for i ∈ 1:length(midplane_cell_centers)] if length(r) != length(quantity) - throw(DimensionMismatch(string( - "Number of cell center coordinates (", length(r), - ") does not match number of cell values (", length(quantity), - ") in the midplane subset." - ))) + throw( + DimensionMismatch( + string( + "Number of cell center coordinates (", length(r), + ") does not match number of cell values (", length(quantity), + ") in the midplane subset.", + ), + ), + ) end eq_idx = 1 - r_eq = dd.equilibrium.time_slice[eq_it].profiles_2d[eq_idx].grid.dim1 - z_eq = dd.equilibrium.time_slice[eq_it].profiles_2d[eq_idx].grid.dim2 - rho1_eq = dd.equilibrium.time_slice[eq_it].profiles_1d.rho_tor_norm - psia = dd.equilibrium.time_slice[eq_it].global_quantities.psi_axis - psib = dd.equilibrium.time_slice[eq_it].global_quantities.psi_boundary - psi1_eq = (dd.equilibrium.time_slice[eq_it].profiles_1d.psi .- psia) ./ (psib - psia) - psi2_eq = (dd.equilibrium.time_slice[eq_it].profiles_2d[eq_idx].psi .- psia) ./ (psib - psia) + r_eq = dd.equilibrium.time_slice[eq_time_idx].profiles_2d[eq_idx].grid.dim1 + z_eq = dd.equilibrium.time_slice[eq_time_idx].profiles_2d[eq_idx].grid.dim2 + rho1_eq = dd.equilibrium.time_slice[eq_time_idx].profiles_1d.rho_tor_norm + psia = dd.equilibrium.time_slice[eq_time_idx].global_quantities.psi_axis + psib = dd.equilibrium.time_slice[eq_time_idx].global_quantities.psi_boundary + psi1_eq = + (dd.equilibrium.time_slice[eq_time_idx].profiles_1d.psi .- psia) ./ + (psib - psia) + psi2_eq = + (dd.equilibrium.time_slice[eq_time_idx].profiles_2d[eq_idx].psi .- psia) ./ + (psib - psia) println(size(psi2_eq), ", ", size(r_eq), ", ", size(z_eq)) rzpi = Interpolations.LinearInterpolation((z_eq, r_eq), psi2_eq) - in_bounds = (r .< maximum(r_eq)) .& (r .> minimum(r_eq)) .& (z .> minimum(z_eq)) .& (z .< maximum(z_eq)) + in_bounds = + (r .< maximum(r_eq)) .& (r .> minimum(r_eq)) .& (z .> minimum(z_eq)) .& + (z .< maximum(z_eq)) println(in_bounds) psi_for_quantity = 10.0 .+ zeros(length(r)) psi_for_quantity[in_bounds] = rzpi.(z[in_bounds], r[in_bounds]) @@ -140,22 +190,183 @@ function fill_in_extrapolated_core_profile(dd::OMAS.dd, quantity_name::String) drho = diff(rho1_eq) prepend!(dpsi, [0.0]) prepend!(drho, [0.0]) - rho_for_quantity[in_bounds] = Interpolations.LinearInterpolation(psi1_eq, rho1_eq).(psi_for_quantity[in_bounds]) - + rho_for_quantity[in_bounds] = + Interpolations.LinearInterpolation( + psi1_eq, + rho1_eq, + ).( + 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. - dd.core_profiles.profiles_1d[it].grid.rho_tor_norm = collect(LinRange(0, 1, 201)) + # 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 rho_core = dd.core_profiles.profiles_1d[it].grid.rho_tor_norm # Finally, we're ready to call the extrapolation function and write the result - quantity_core = extrapolate_core(rho_for_quantity, quantity, rho_core) + if method == "simple" + quantity_core = extrapolate_core(rho_for_quantity, quantity, rho_core) + else + throw(ArgumentError(string( + "Unrecognized extraplation method: ", method, + ))) + end parent = dd.core_profiles.profiles_1d[it] tags = split(quantity_name, ".") - for tag in tags[1:end-1] + for tag ∈ tags[1:end-1] parent = getproperty(parent, Symbol(tag)) end - setproperty!(parent, Symbol(tags[end]), quantity_core) -end \ No newline at end of file + return setproperty!(parent, Symbol(tags[end]), quantity_core) +end + +""" + function extrapolate_edge_exp( + quantity_edge::Vector{Float64}, + dqdpsi::Vector{Float64}, + psin_out::Vector{Float64}, + ) + +Exponential decay version of edge profile extrapolation. Should work well for +many quantities, including Te and ne. If the exponential profile has no background +offset, then its amplitude and scale length can be completely defined by matching +the quantity value and first derivative at the edge of the mesh. + +quantity_edge: values of some physics quantity in cells along the outer edge of the mesh +dqdpsi: Gradient of the quantity vs. psi, aligned perpendicular to the row of cells +being used. +psin_out: Normalized psi values along a vector orthogonal to the row of cells along the +edge. These psi_N values should be outside of the mesh (because the quantity is +already known in the mesh). +The output will be a matrix. +""" +function extrapolate_edge_exp( + quantity_edge::Vector{Float64}, + dqdpsi::Vector{Float64}, + psin_out::Vector{Float64}, +) + x = psin_out - 1.0 + lambda = -quantity_edge / dqdpsi + q0 = quantity_edge ./ exp(-x' ./ lambda) + return y0 * exp(-x ./ lambda) +end + +""" + function mesh_psi_spacing(dd::OMAS.dd; eq_time_idx::Int64=1) + +Inspects the mesh to see how far apart faces are in psi_N. +Requires that GGD and equilibrium are populated. + +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. +""" +function mesh_psi_spacing(dd::OMAS.dd; eq_time_idx::Int64=1) + # Inspect input + if length(dd.equilibrium.time_slice) < eq_time_idx + throw( + ArgumentError( + string( + "DD equilibrium does not have enough time slices: ", + length(dd.equilibrium.time_slice), + " slices were present but slice index ", + eq_time_idx, " was requested.", + ), + ), + ) + end + bad_ggd = (length(dd.edge_profiles.grid_ggd) < 1) + if bad_ggd + throw(ArgumentError(string( + "Invalid GGD data.", + ))) + end + + # Get flux map + eq_idx = 1 # Most cases won't have need for more than 1 of these, as far as I know + eqt = dd.equilibrium.time_slice[eq_time_idx] + p2 = eqt.profiles_2d[eq_idx] + r_eq = p2.grid.dim1 + z_eq = p2.grid.dim2 + psi = p2.psi + psia = eqt.global_quantities.psi_axis + psib = eqt.global_quantities.psi_boundary + psin_eq = (psi .- psia) ./ (psib - psia) + rzpi = Interpolations.LinearInterpolation((z_eq, r_eq), psin_eq) + println(minimum(r_eq), ", ", maximum(r_eq)) + println(minimum(z_eq), ", ", maximum(z_eq)) + + # Get a row of cells. Since the mesh should be aligned to the flux surfaces, + # it shouldn't matter which row is used, although the divertor rows might be + # weird. So use the outboard midplane. That's always a solid choice. + ggd_idx = 1 + space_number = 1 + space = dd.edge_profiles.grid_ggd[ggd_idx].space[space_number] + midplane_subset = + SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[ggd_idx], 11) + midplane_cell_centers = GGDUtils.get_subset_centers(space, midplane_subset) + r_mesh = [midplane_cell_centers[i][1] for i ∈ 1:length(midplane_cell_centers)] + z_mesh = [midplane_cell_centers[i][2] for i ∈ 1:length(midplane_cell_centers)] + println(minimum(r_mesh), ", ", maximum(r_mesh)) + println(minimum(z_mesh), ", ", maximum(z_mesh)) + psin_mesh = rzpi.(z_mesh, r_mesh) + # This should come out sorted, but with GGD, who knows. + ii = sortperm(psin_mesh) + psin = psin_mesh[ii] + dpsin = psin[end] - psin[end-1] + return dpsin +end + +""" + fill_in_extrapolated_edge_profile!( + dd::OMAS.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, + quantity_name::String; + method::String="simple", + eq_time_idx::Int64=1, +) + # Inspect input + if length(dd.equilibrium.time_slice) < eq_time_idx + throw( + ArgumentError( + string( + "DD equilibrium does not have enough time slices: ", + length(dd.equilibrium.time_slice), + " slices were present but slice index ", + eq_time_idx, " was requested.", + ), + ), + ) + end + bad_ggd = (length(dd.edge_profiles.grid_ggd) < 1) + if bad_ggd + throw(ArgumentError(string( + "Invalid GGD data.", + ))) + end + + if method == "simple" + println("THERE IS NOT ACTUALLY AN EDGE EXTRAPOLATION METHOD YET.") + println("THIS IS A PLACEHOLDER SO THE REST OF THE WORKFLOW CAN BE SET UP.") + # In this case, there's a good chance that we'll need a few different classes + # of extrapolation techniques for different quantities as the may behave + # differently. + else + throw(ArgumentError(string( + "Unrecognized extraplation method: ", method, + ))) + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 943b396..10644b1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,7 +1,7 @@ -import SD4SOLPS -import SOLPS2IMAS -import OMAS -import EFIT +using SD4SOLPS: SD4SOLPS +using SOLPS2IMAS: SOLPS2IMAS +using OMAS: OMAS +using EFIT: EFIT using Plots using Test @@ -12,20 +12,29 @@ Makes a modified tanh profile based on the formula from: Groebner, et al., Nucl. Fusion 41, 1789 (2001) 10.1088/0029-5515/41/12/306 x: position basis. The output will be a function of this coordinate. - Could be psi_N or R or Z or whatever, as long as sym and hwid are - specified in terms of the same coordinate. - Defaults are given in terms of psi_N. +Could be psi_N or R or Z or whatever, as long as sym and hwid are +specified in terms of the same coordinate. +Defaults are given in terms of psi_N. sym: location of the tanh symmetry point hwid: half width of the tanh offset: value at the foot of the tanh / SOL value pedestal: value at the top of the tanh / pedestal value alpha: interior slope """ -function make_test_profile(x; sym=0.9657, hwid=0.05421, offset=0.0953, pedestal=3.58901, alpha=0.005) +function make_test_profile( + x; + sym=0.9657, + hwid=0.05421, + offset=0.0953, + pedestal=3.58901, + alpha=0.005, +) z = (sym .- x) ./ hwid amplitude = (pedestal - offset) / 2.0 b = offset + amplitude - return b .+ amplitude .* ((1 .+ alpha .* z) .* exp.(z) .- exp.(.-z)) ./ (exp.(z) .+ exp.(.-z)) + return b .+ + amplitude .* ((1 .+ alpha .* z) .* exp.(z) .- exp.(.-z)) ./ + (exp.(z) .+ exp.(.-z)) end # function plot_test_profiles() @@ -37,6 +46,30 @@ end # Plots.plot!(edge_rho, edge_quantity, marker='.') # end +function define_default_sample_set() + sample_path = + splitdir(pathof(SD4SOLPS))[1] * + "/../sample/ITER_Lore_2296_00000/extended_output" + sample_path2 = + splitdir(pathof(SD4SOLPS))[1] * "/../sample/ITER_Lore_2296_00000/baserun" + sample_path3 = + splitdir(pathof(SD4SOLPS))[1] * "/../sample/ITER_Lore_2296_00000/run_restart" + sample_path4 = splitdir(pathof(SOLPS2IMAS))[1] * "/../samples" + + # Requires dvc pull before the full samples will be loaded + # Sorry, the minimal samples are too minimal for this one. + file_list = SD4SOLPS.find_files_in_allowed_folders( + sample_path, sample_path2, sample_path3, sample_path4; + eqdsk_file="thereisntoneyet", + allow_reduced_versions=false, + ) + b2fgmtry, b2time, b2mn, gridspec, eqdsk = file_list + eqdsk = + splitdir(pathof(SD4SOLPS))[1] * + "/../sample/ITER_Lore_2296_00000/EQDSK/Baseline2008-li0.70.x4.mod2.eqdsk" + return b2fgmtry, b2time, b2mn, gridspec, eqdsk +end + @testset "core_profile_extension" begin # Just the basic profile extrapolator ------------------ edge_rho = Array(LinRange(0.88, 1.0, 18)) @@ -47,20 +80,7 @@ end # The full workflow -------------------------------------- # Setup sample DD - sample_path = splitdir(pathof(SD4SOLPS))[1] * "/../sample/ITER_Lore_2296_00000/extended_output" - sample_path2 = splitdir(pathof(SD4SOLPS))[1] * "/../sample/ITER_Lore_2296_00000/baserun" - sample_path3 = splitdir(pathof(SD4SOLPS))[1] * "/../sample/ITER_Lore_2296_00000/run_restart" - sample_path4 = splitdir(pathof(SOLPS2IMAS))[1] * "/../samples" - - # Requires dvc pull before the full samples will be loaded - # Sorry, the minimal samples are too minimal for this one. - file_list = SD4SOLPS.find_files_in_allowed_folders( - sample_path, sample_path2, sample_path3, sample_path4, - eqdsk_file="thereisntoneyet", - allow_reduced_versions=false, - ) - b2fgmtry, b2time, b2mn, gridspec, eqdsk = file_list - eqdsk = splitdir(pathof(SD4SOLPS))[1] * "/../sample/ITER_Lore_2296_00000/EQDSK/Baseline2008-li0.70.x4.mod2.eqdsk" + b2fgmtry, b2time, b2mn, gridspec, eqdsk = define_default_sample_set() println(b2fgmtry) println(b2time) println(b2mn) @@ -72,11 +92,11 @@ end @test isfile(b2mn) @test isfile(gridspec) @test isfile(eqdsk) - dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) + dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) SD4SOLPS.geqdsk_to_imas(eqdsk, dd) rho = dd.equilibrium.time_slice[1].profiles_1d.rho_tor_norm - - if !SD4SOLPS.check_rho_1d(dd, time_slice=1) + + if !SD4SOLPS.check_rho_1d(dd; time_slice=1) SD4SOLPS.add_rho_to_equilibrium!(dd) rho = dd.equilibrium.time_slice[1].profiles_1d.rho_tor_norm println("Repaired missing rho for core profile test") @@ -88,27 +108,43 @@ end test_slice_idx = 1 # Do it - SD4SOLPS.fill_in_extrapolated_core_profile(dd, quantity_name) + SD4SOLPS.fill_in_extrapolated_core_profile!(dd, quantity_name) # Inspect results @test length(dd.core_profiles.profiles_1d) > 0 core_prof = dd.core_profiles.profiles_1d[test_slice_idx] tags = split(quantity_name, ".") quantity = dd.core_profiles.profiles_1d[test_slice_idx] - for tag in tags + for tag ∈ tags quantity = getproperty(quantity, Symbol(tag)) end rho_core = dd.core_profiles.profiles_1d[test_slice_idx].grid.rho_tor_norm @test length(quantity) > 0 @test length(quantity) == length(rho_core) end + +@testset "edge_profile_extension" begin + # Test for getting mesh spacing + b2fgmtry, b2time, b2mn, gridspec, eqdsk = define_default_sample_set() + dd = SOLPS2IMAS.solps2imas(b2fgmtry, b2time, gridspec, b2mn) + SD4SOLPS.geqdsk_to_imas(eqdsk, dd) + dpsin = SD4SOLPS.mesh_psi_spacing(dd) + @test dpsin > 0.0 + + n_edge = 47 + n_outer_prof = 13 + quantity_edge = [collect(LinRange(2.0, 75.0, 7)); collect(LinRange(75.0, 98.3, 17))] + quantity_edge = [quantity_edge; reverse(quantity_edge)[2:end]] + gradient_edge = [collect(LinRange(-1, -100, 7)); collect(LinRange(-100, -350, 17))] + gradient_edge = [gradient_edge; reverse(gradient_edge)[2:end]] + psin_out = collect(LinRange(1.0, 1.25, n_outer_prof + 1))[2:end] end @testset "utilities" begin # Test for finding files in allowed folders sample_path = splitdir(pathof(SOLPS2IMAS))[1] * "/../samples/" file_list = SD4SOLPS.find_files_in_allowed_folders( - sample_path, eqdsk_file="thereisntoneyet", allow_reduced_versions=true, + sample_path; eqdsk_file="thereisntoneyet", allow_reduced_versions=true, ) @test length(file_list) == 5 b2fgmtry, b2time, b2mn, gridspec, eqdsk = file_list @@ -121,7 +157,8 @@ end @test length(gridspec) > 10 @test endswith(gridspec, "gridspacedesc.yml") - # Test for sweeping 1D core profiles into 2D R,Z (or anyway evaluating them at any R,Z location) + # Test for sweeping 1D core profiles into 2D R,Z + # (or anyway evaluating them at any R,Z location) dd = OMAS.dd() eqdsk_file = splitdir(pathof(SD4SOLPS))[1] * "/../sample/geqdsk_iter_small_sample" SD4SOLPS.geqdsk_to_imas(eqdsk_file, dd) @@ -133,7 +170,8 @@ end 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) + dd.core_profiles.profiles_1d[prof_time_idx].electrons.density = + make_test_profile(rho_n) r_mag_axis = dd.equilibrium.time_slice[1].global_quantities.magnetic_axis.r z_mag_axis = dd.equilibrium.time_slice[1].global_quantities.magnetic_axis.z rg = r_mag_axis .* [0.85, 0.90, 0.95, 1.0, 1.05, 1.10, 1.15] @@ -141,11 +179,12 @@ end points = collect(Base.Iterators.product(rg, zg)) r = getfield.(points, 1) z = getfield.(points, 2) - if !SD4SOLPS.check_rho_1d(dd, time_slice=eq_time_idx) + if !SD4SOLPS.check_rho_1d(dd; time_slice=eq_time_idx) SD4SOLPS.add_rho_to_equilibrium!(dd) println("DD was repaired (rho added) for core 2d utility test") end - density_on_grid = SD4SOLPS.core_profile_2d(dd, prof_time_idx, eq_time_idx, quantity, r, z) + density_on_grid = + SD4SOLPS.core_profile_2d(dd, prof_time_idx, eq_time_idx, quantity, r, z) @test size(density_on_grid) == (length(rg), length(zg)) end @@ -156,7 +195,7 @@ end SD4SOLPS.geqdsk_to_imas(eqdsk, dd) # Make sure rho is missing nt = length(dd.equilibrium.time_slice) - for it = 1:nt + for it ∈ 1:nt eqt = dd.equilibrium.time_slice[it] eqt.profiles_1d.rho_tor_norm = Vector{Float64}() end @@ -169,13 +208,14 @@ end end @testset "geqdsk_to_imas" begin - sample_files = (splitdir(pathof(SD4SOLPS))[1] * "/../sample/") .* [ - "g184833.03600", "geqdsk_iter_small_sample" - ] + sample_files = + (splitdir(pathof(SD4SOLPS))[1] * "/../sample/") .* [ + "g184833.03600", "geqdsk_iter_small_sample", + ] tslice = 1 - for sample_file in sample_files + for sample_file ∈ sample_files dd = OMAS.dd() - SD4SOLPS.geqdsk_to_imas(sample_file, dd, time_index=tslice) + SD4SOLPS.geqdsk_to_imas(sample_file, dd; time_index=tslice) eqt = dd.equilibrium.time_slice[tslice] # global @@ -215,13 +255,27 @@ end @testset "preparation" begin eqdsk_file = "geqdsk_iter_small_sample" sample_paths = [ - splitdir(pathof(SOLPS2IMAS))[1] * "/../samples/", splitdir(pathof(SD4SOLPS))[1] * "/../sample/", + splitdir(pathof(SOLPS2IMAS))[1] * "/../samples/", ] - dd = SD4SOLPS.preparation(eqdsk_file, sample_paths...) + core_method = "simple" + edge_method = "simple" + filename = splitdir(pathof(SD4SOLPS))[1] * "/../sd_input_data" + output_format = "json" + dd = SD4SOLPS.preparation( + eqdsk_file, + sample_paths...; + core_method=core_method, + edge_method=edge_method, + filename=filename, + output_format=output_format, + ) + out_file = filename * "." * output_format p2 = dd.equilibrium.time_slice[1].profiles_2d[1] psirz = p2.psi r = p2.grid.dim1 z = p2.grid.dim2 @test size(psirz) == (length(r), length(z)) + println(out_file) + @test isfile(out_file) end