From 7f73bc7f3e4752dcd1b6984dbc64b76c131048e4 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Thu, 21 Sep 2023 12:28:35 -0700 Subject: [PATCH] Formatted with config file and adding format check * Added format file .JuliaFormatter * Added GitHub Action for format checking * Formatted the code with the file above and some hand tweaks on comments --- .JuliaFormatter.toml | 27 ++ .github/workflows/format_check.yml | 36 ++ src/SOLPS2IMAS.jl | 717 ++++++++++++++++++----------- test/runtests.jl | 56 ++- 4 files changed, 552 insertions(+), 284 deletions(-) create mode 100644 .JuliaFormatter.toml create mode 100644 .github/workflows/format_check.yml 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/src/SOLPS2IMAS.jl b/src/SOLPS2IMAS.jl index 2e0f3a0..130001b 100644 --- a/src/SOLPS2IMAS.jl +++ b/src/SOLPS2IMAS.jl @@ -1,7 +1,7 @@ module SOLPS2IMAS using Revise -import OMAS +using OMAS: OMAS using NCDatasets: Dataset, dimnames using YAML: load_file as YAML_load_file export try_omas @@ -10,7 +10,6 @@ export read_b2_output export search_points export solps2imas - function try_omas() ids = OMAS.dd() resize!(ids.equilibrium.time_slice, 1) @@ -18,7 +17,6 @@ function try_omas() return nothing end - function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69) """This doesn't work well""" center_r = 1.7 @@ -27,9 +25,9 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69) inner_minor = 0.3 outer_minor = 0.8 elongation = 1.8 - θₓ = 8 * π/6 - θᵢ = θₓ - π/4 - θₒ = θₓ + π/4 + θₓ = 8 * π / 6 + θᵢ = θₓ - π / 4 + θₒ = θₓ + π / 4 inner_leg_length = 0.5 outer_leg_length = 0.6 dx_il = inner_leg_length / lcut @@ -42,8 +40,8 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69) cry = Array{Float64}(undef, 4, ny, nx) xpoint_r = center_r + aminor * cos(θₓ) xpoint_z = center_z + aminor * sin(θₓ) * elongation - for ix::Int = 1:nx - for iy::Int = 1:ny + for ix::Int ∈ 1:nx + for iy::Int ∈ 1:ny if iy < sep dy = dr_in else @@ -57,11 +55,15 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69) cell_centers_r[iy, ix] = cos(θ) * a + center_r cell_centers_z[iy, ix] = sin(θ) * a * elongation + center_z elseif ix < lcut - cell_centers_r[iy, ix] = xpoint_r + cos(θᵢ) * dx_il * (lcut - ix) + sin(θᵢ) * dy * (iy - sep) - cell_centers_z[iy, ix] = xpoint_z + sin(θᵢ) * dx_il * (lcut - ix) + cos(θᵢ) * dy * (iy - sep) + cell_centers_r[iy, ix] = + xpoint_r + cos(θᵢ) * dx_il * (lcut - ix) + sin(θᵢ) * dy * (iy - sep) + cell_centers_z[iy, ix] = + xpoint_z + sin(θᵢ) * dx_il * (lcut - ix) + cos(θᵢ) * dy * (iy - sep) else - cell_centers_r[iy, ix] = xpoint_r + cos(θₒ) * dx_ol * (ix - rcut) + sin(θₒ) * dy * (iy - sep) - cell_centers_z[iy, ix] = xpoint_z + sin(θₒ) * dx_ol * (ix - rcut) + cos(θₒ) * dy * (iy - sep) + cell_centers_r[iy, ix] = + xpoint_r + cos(θₒ) * dx_ol * (ix - rcut) + sin(θₒ) * dy * (iy - sep) + cell_centers_z[iy, ix] = + xpoint_z + sin(θₒ) * dx_ol * (ix - rcut) + cos(θₒ) * dy * (iy - sep) end end end @@ -70,7 +72,7 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69) rmax = 3.0 zmin = -3.0 zmax = 3.0 - for i::Int = 1:nx*ny + for i::Int ∈ 1:nx*ny # if mod(i, nx) == 0 # print(cell_centers_r[i], " ", cell_centers_r[i] > rmax, " ", rmax) # end @@ -85,7 +87,6 @@ function generate_test_data(nx=94, ny=38, sep=15, lcut=25, rcut=69) return (cell_centers_r, cell_centers_z) end - function read_b2time_output(filename) dim_order = ( "time", @@ -98,13 +99,17 @@ function read_b2time_output(filename) ) ret_dict = Dict("dim" => Dict(), "data" => Dict()) ds = Dataset(filename) - for key in keys(ds.dim) + for key ∈ keys(ds.dim) ret_dict["dim"][key] = ds.dim[key] end - for key in keys(ds) + for key ∈ keys(ds) if key != "ntstep" d = dimnames(ds[key]) - permute = [y for y in [findfirst(x->x==dimord, d) for dimord in dim_order] if y !== nothing] + permute = [ + y for + y ∈ [findfirst(x -> x == dimord, d) for dimord ∈ dim_order] if + y !== nothing + ] try ret_dict["data"][key] = permutedims(Array(ds[key]), permute) catch e @@ -117,13 +122,12 @@ function read_b2time_output(filename) return ret_dict end - function read_b2mn_output(filename) lines = open(filename) do f - readlines(f) + return readlines(f) end contents = Dict() - for line in lines + for line ∈ lines if startswith(line, "'") splits = split(line, "'"; keepempty=false) contents[splits[1]] = parse(Float64, splits[end]) @@ -132,10 +136,7 @@ function read_b2mn_output(filename) return contents end - function read_b2_output(filename) - # x = readdlm(filename, ' ', Float64, '\n', header=true) # Doesn't handle the text lines at the start of each array - if cmp(splitext(filename)[2], ".nc") == 0 return read_b2time_output(filename) end @@ -144,7 +145,7 @@ function read_b2_output(filename) array_sizes = Dict() ret_dict = Dict() lines = open(filename) do f - readlines(f) + return readlines(f) end nx = 0 ny = 0 @@ -153,7 +154,7 @@ function read_b2_output(filename) arraysize = 0 arraytype = nothing j = 1 - for l in lines + for l ∈ lines if startswith(l, "*cf:") j = 1 # Reset intra-array element counter _, arraytype, arraysize, tag = split(l) @@ -168,10 +169,10 @@ function read_b2_output(filename) array_sizes[tag] = arraysize elseif tag != "" if arraytype == "int" - array_line = [parse(Int, ss) for ss in split(l)] + array_line = [parse(Int, ss) for ss ∈ split(l)] array_inc = size(array_line)[1] elseif arraytype == "real" - array_line = [parse(Float64, ss) for ss in split(l)] + array_line = [parse(Float64, ss) for ss ∈ split(l)] array_inc = size(array_line)[1] else array_line = l @@ -190,12 +191,14 @@ function read_b2_output(filename) elseif "nx,ny,ns" ∈ keys(contents) return extract_state_quantities(contents) else - throw(DomainError(keys(contents), - "nx,ny (b2fgmtry) or nx,ny,ns (b2fstate) must be present in b2 output file.")) + throw( + DomainError(keys(contents), + "nx,ny (b2fgmtry) or nx,ny,ns (b2fstate) must be present in b2 output file.", + ), + ) end end - function extract_geometry(gmtry) ret_dict = Dict("dim" => Dict(), "data" => Dict()) ret_dict["dim"]["nx_no_guard"], ret_dict["dim"]["ny_no_guard"] = gmtry["nx,ny"] @@ -208,12 +211,13 @@ function extract_geometry(gmtry) # Adding placeholder timestamp ret_dict["dim"]["time"] = 1 ret_dict["data"]["timesa"] = [0.0] - for k in keys(gmtry) + for k ∈ keys(gmtry) # The 4 fields of bb are poloidal, radial, toroidal, and total magnetic field # according to page 212 of D. Coster, "SOLPS-ITER [manual]" (2019) # The 4 fields in crx and cry are the corners of each grid cell. if k ∈ ["crx", "cry", "bb"] - ret_dict["data"][k] = permutedims(reshape(gmtry[k], (nx, ny, 4, 1)), (4, 3, 2, 1)) + ret_dict["data"][k] = + permutedims(reshape(gmtry[k], (nx, ny, 4, 1)), (4, 3, 2, 1)) elseif k ∈ ["leftcut", "bottomcut", "rightcut", "topcut"] ret_dict["data"][k] = Array([gmtry[k][1]]) elseif k ∈ ["leftcut2", "bottomcut2", "rightcut2", "topcut2"] @@ -227,10 +231,11 @@ function extract_geometry(gmtry) return ret_dict end - function extract_state_quantities(state) ret_dict = Dict("dim" => Dict(), "data" => Dict()) - ret_dict["dim"]["nx_no_guard"], ret_dict["dim"]["ny_no_guard"], ret_dict["dim"]["ns"] = state["nx,ny,ns"] + ret_dict["dim"]["nx_no_guard"], + ret_dict["dim"]["ny_no_guard"], + ret_dict["dim"]["ns"] = state["nx,ny,ns"] # includes guard cells nx = ret_dict["dim"]["nx"] = ret_dict["dim"]["nx_no_guard"] + 2 ny = ret_dict["dim"]["ny"] = ret_dict["dim"]["ny_no_guard"] + 2 @@ -239,16 +244,19 @@ function extract_state_quantities(state) # Adding placeholder timestamp ret_dict["dim"]["time"] = 1 ret_dict["data"]["timesa"] = [0.0] - for k in keys(state) + for k ∈ keys(state) l = length(state[k]) if l == nx * ny ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, 1)), (3, 2, 1)) elseif l == nx * ny * ns - ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, ns, 1)), (4, 3, 2, 1)) + ret_dict["data"][k] = + permutedims(reshape(state[k], (nx, ny, ns, 1)), (4, 3, 2, 1)) elseif l == nx * ny * ndir - ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, ndir, 1)), (4, 3, 2, 1)) + ret_dict["data"][k] = + permutedims(reshape(state[k], (nx, ny, ndir, 1)), (4, 3, 2, 1)) elseif l == nx * ny * ndir * ns - ret_dict["data"][k] = permutedims(reshape(state[k], (nx, ny, ndir, ns, 1)), (5, 4, 3, 2, 1)) + ret_dict["data"][k] = + permutedims(reshape(state[k], (nx, ny, ndir, ns, 1)), (5, 4, 3, 2, 1)) elseif l == ns ret_dict["data"][k] = permutedims(reshape(state[k], (ns, 1)), (2, 1)) elseif k ∉ keys(ret_dict["dim"]) @@ -258,7 +266,6 @@ function extract_state_quantities(state) return ret_dict end - """ xytoc(ix, iy; nx) @@ -270,7 +277,6 @@ function xytoc(ix, iy; nx) return ic end - """ ctoxy(ic; nx) @@ -282,9 +288,8 @@ function ctoxy(ic; nx) return ix, iy end - """ - in_core(; ix, iy, topcut, bottomcut, leftcut, rightcut) +in_core(; ix, iy, topcut, bottomcut, leftcut, rightcut) Returns true if cell indexed ix, iy lie inside the core """ @@ -292,16 +297,12 @@ function in_core(; ix, iy, topcut, bottomcut, leftcut, rightcut) return bottomcut + 1 < iy < topcut + 2 && leftcut + 1 < ix < rightcut + 2 end - """ in_sol(; iy, topcut, kwargs...) Returns true if cell indexed ix, iy lie inside the SOL """ -function in_sol(; iy, topcut, kwargs...) - return topcut + 1 < iy -end - +in_sol(; iy, topcut, kwargs...) = topcut + 1 < iy """ in_idr(; ix, iy, topcut, bottomcut, leftcut, kwargs...) @@ -312,7 +313,6 @@ function in_idr(; ix, iy, topcut, bottomcut, leftcut, kwargs...) return bottomcut + 1 < iy < topcut + 2 && ix < leftcut + 2 end - """ in_odr(; ix, iy, topcut, bottomcut, rightcut, kwargs...) @@ -328,52 +328,66 @@ end # Thus, x-edges will have odd boundary index and y_edges will have even # List of tuples (boundary_ind, (corner pair forming edge)) chosen_edge_order = [(1, (1, 2)), - (2, (2, 4)), - (3, (4, 3)), - (4, (3, 1))] - + (2, (2, 4)), + (3, (4, 3)), + (4, (3, 1))] """ is_x_aligned(;boundary_ind) y_aligned edges will have odd boundary_ind based on chosen order of numbering them """ -function is_x_aligned(;boundary_ind) - return mod(boundary_ind, 2) == 1 -end - +is_x_aligned(; boundary_ind) = mod(boundary_ind, 2) == 1 """ is_y_aligned(; boundary_ind) y_aligned edges will have even boundary_ind based on chosen order of numbering them """ -function is_y_aligned(; boundary_ind) - return mod(boundary_ind, 2) == 0 -end +is_y_aligned(; boundary_ind) = mod(boundary_ind, 2) == 0 """ - is_core_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut) +is_core_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut) Returns true if boundary_ind of a cell at ix, iy is on core_cut (Y-aliged edge) """ -function is_core_cut(; ix, iy, cells, nx, boundary_ind, topcut, bottomcut, leftcut, rightcut) +function is_core_cut(; + ix, + iy, + cells, + nx, + boundary_ind, + topcut, + bottomcut, + leftcut, + rightcut, +) if bottomcut + 1 < iy < topcut + 2 && ix == leftcut + 2 && mod(boundary_ind, 2) == 0 ixr = rightcut + 1 this_cell = cells[xytoc(ix, iy; nx=nx)] - other_side_cell_ind = xytoc(ixr, iy; nx=nx) # Cell on the other side of core cut + # Cell on the other side of core cut + other_side_cell_ind = xytoc(ixr, iy; nx=nx) return other_side_cell_ind ∈ this_cell.boundary[boundary_ind].neighbours end return false end - """ - is_pfr_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut) +is_pfr_cut(; ix, iy, nx, cells, boundary_ind, topcut, bottomcut, leftcut, rightcut) Returns true if boundary_ind of a cell at ix, iy is on core_cut (Y-aliged edge) """ -function is_pfr_cut(; ix, iy, cells, nx, boundary_ind, topcut, bottomcut, leftcut, rightcut) +function is_pfr_cut(; + ix, + iy, + cells, + nx, + boundary_ind, + topcut, + bottomcut, + leftcut, + rightcut, +) if bottomcut + 1 < iy < topcut + 2 && ix == leftcut + 1 && mod(boundary_ind, 2) == 0 ixr = rightcut + 2 this_cell = cells[xytoc(ix, iy; nx=nx)] @@ -383,7 +397,6 @@ function is_pfr_cut(; ix, iy, cells, nx, boundary_ind, topcut, bottomcut, leftc return false end - """ is_outer_throat(; ix, iy, boundary_ind, topcut, rightcut, kwargs...) @@ -393,7 +406,6 @@ function is_outer_throat(; ix, iy, boundary_ind, topcut, rightcut, kwargs...) return topcut + 1 < iy && ix == rightcut + 1 && boundary_ind == 2 end - """ is_inner_throat(; ix, iy, boundary_ind, topcut, leftcut, kwargs...) @@ -403,7 +415,6 @@ function is_inner_throat(; ix, iy, boundary_ind, topcut, leftcut, kwargs...) return topcut + 1 < iy && ix == leftcut + 2 && boundary_ind == 4 end - """ is_outer_midplane(; ix, jxa, boundary_ind) @@ -414,7 +425,6 @@ function is_outer_midplane(; ix, iy, jxa, boundary_ind, topcut, kwargs...) return ix == jxa && boundary_ind == 2 end - """ is_inner_midplane(; ix, jxa, boundary_ind) @@ -425,37 +435,47 @@ function is_inner_midplane(; ix, iy, jxi, boundary_ind, topcut, kwargs...) return ix == jxi && boundary_ind == 4 end - """ - is_outer_target(; ix, nx, boundary_ind, kwargs...) + is_outer_target(; ix, nx, boundary_ind) Returns true if boundary_ind of a cell at ix, iy is on outer target """ -function is_outer_target(; ix, nx, boundary_ind) - return ix == nx && boundary_ind == 2 -end - +is_outer_target(; ix, nx, boundary_ind) = ix == nx && boundary_ind == 2 """ is_inner_target(; ix, boundary_ind, kwargs...) Returns true if boundary_ind of a cell at ix, iy is on inner target """ -function is_inner_target(; ix, boundary_ind) - return ix == 1 && boundary_ind == 4 -end - +is_inner_target(; ix, boundary_ind) = ix == 1 && boundary_ind == 4 """ - is_core_boundary(; ix, iy, boundary_ind, bottomcut, leftcut, rightcut, kwargs...) - -Returns true if boundary_ind of a cell at ix, iy is on core boundary (central blank spot boundary) + is_core_boundary(; + ix, + iy, + boundary_ind, + bottomcut, + leftcut, + rightcut, + kwargs..., + +) + +Returns true if boundary_ind of a cell at ix, iy is on core boundary (central blank +spot boundary) """ -function is_core_boundary(; ix, iy, boundary_ind, bottomcut, leftcut, rightcut, kwargs...) +function is_core_boundary(; + ix, + iy, + boundary_ind, + bottomcut, + leftcut, + rightcut, + kwargs..., +) return bottomcut + 2 == iy && leftcut + 1 < ix < rightcut + 2 && boundary_ind == 1 end - """ is_separatix(; iy, boundary_ind, topcut, kwargs...) @@ -465,7 +485,6 @@ function is_separatix(; iy, boundary_ind, topcut, kwargs...) return topcut + 2 == iy && boundary_ind == 1 end - """ add_subset_element!(subset, sn, dim, index::Int, in_subset=(x...)->true; kwargs...) @@ -474,7 +493,14 @@ with space number sn, dimension dim, index. To determine, if the element should be added or not, a function in_subset can be provided that gets the arguments (kwargs...). These functions will be in_core, in_sol etc as difined above. """ -function add_subset_element!(subset, sn, dim, index::Int, in_subset=(x...)->true; kwargs...) +function add_subset_element!( + subset, + sn, + dim, + index::Int, + in_subset=(x...) -> true; + kwargs..., +) if in_subset(; kwargs...) dd_ind = length(subset.element) + 1 resize!(subset.element, dd_ind) @@ -485,18 +511,32 @@ function add_subset_element!(subset, sn, dim, index::Int, in_subset=(x...)->true end end - """ - add_subset_element!(subset, sn, dim, index::Vector{Int}, in_subset=(x...)->true; kwargs...) + add_subset_element!( + subset, + sn, + dim, + index::Vector{Int}, + in_subset=(x...) -> true; + kwargs..., + +) Overloaded to work differently (faster) with list of indices to be added. """ -function add_subset_element!(subset, sn, dim, index::Vector{Int}, in_subset=(x...)->true; kwargs...) +function add_subset_element!( + subset, + sn, + dim, + index::Vector{Int}, + in_subset=(x...) -> true; + kwargs..., +) if in_subset(; kwargs...) dd_start_ind = length(subset.element) + 1 resize!(subset.element, length(subset.element) + length(index)) dd_stop_ind = length(subset.element) - for (ii, dd_ind) in enumerate(dd_start_ind:dd_stop_ind) + for (ii, dd_ind) ∈ enumerate(dd_start_ind:dd_stop_ind) resize!(subset.element[dd_ind].object, 1) subset.element[dd_ind].object[1].space = sn subset.element[dd_ind].object[1].dimension = dim @@ -505,7 +545,6 @@ function add_subset_element!(subset, sn, dim, index::Vector{Int}, in_subset=(x.. end end - """ get_xpoint_nodes(gmtry) @@ -517,23 +556,27 @@ function get_xpoint_nodes(gmtry) nt = gmtry["dim"]["time"] # Find cells around x-point xpcells = [(gmtry["data"]["topcut"][1] + 1, gmtry["data"]["leftcut"][1] + 1), - (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["leftcut"][1] + 2), - (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 1), - (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 2), - (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 1), - (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 2), - (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 1), - (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 2)] + (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["leftcut"][1] + 2), + (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 1), + (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["leftcut"][1] + 2), + (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 1), + (gmtry["data"]["topcut"][1] + 1, gmtry["data"]["rightcut"][1] + 2), + (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 1), + (gmtry["data"]["topcut"][1] + 2, gmtry["data"]["rightcut"][1] + 2)] # Get list of all nodes in these cells candidate_nodes = [] resize!(candidate_nodes, nt) - for it in 1:nt - candidate_nodes[it] = [[[crx[it, icorner, iy, ix], cry[it, icorner, iy, ix]] for icorner in 1:4] for (iy, ix) in xpcells] + for it ∈ 1:nt + candidate_nodes[it] = [ + [ + [crx[it, icorner, iy, ix], cry[it, icorner, iy, ix]] for icorner ∈ 1:4 + ] for (iy, ix) ∈ xpcells + ] end xpoint_nodes = [] resize!(xpoint_nodes, nt) # Find the node that is common among all the cells - for it in 1:nt + for it ∈ 1:nt xpoint_nodes[it] = intersect(candidate_nodes[it]...) end return xpoint_nodes @@ -546,9 +589,10 @@ function search_points(ids, r, z) space_number = 1 subset_idx_node = 1 # If an index remains at 0, it means the point in question was not found - nodes = ids.edge_profiles.grid_ggd[grid_number].space[space_number].objects_per_dimension[subset_idx_node].object - for j = 1:n - for i in eachindex(nodes) + nodes = + ids.edge_profiles.grid_ggd[grid_number].space[space_number].objects_per_dimension[subset_idx_node].object + for j ∈ 1:n + for i ∈ eachindex(nodes) rn = nodes[i].geometry[1] zn = nodes[i].geometry[2] if (rn == r[j]) && (zn == z[j]) @@ -560,36 +604,34 @@ function search_points(ids, r, z) return indices end - """ search_edges(edges, edge_nodes) search if an edge with nodes as edge_nodes already exists """ function search_edges(edges, edge_nodes) - for ii in eachindex(edges) + for ii ∈ eachindex(edges) if edge_nodes[2] == edges[ii].nodes[1] && edge_nodes[1] == edges[ii].nodes[2] return ii - elseif edge_nodes[2] == edges[ii].nodes[1] && edge_nodes[1] == edges[ii].nodes[2] + elseif edge_nodes[2] == edges[ii].nodes[1] && + edge_nodes[1] == edges[ii].nodes[2] return ii end end return 0 end - """ distance_between_nodes(nodes, node_inds) Return distance between two node indices """ function distance_between_nodes(nodes, node_inds) - return √(sum((nodes[node_inds[1]].geometry - nodes[node_inds[2]].geometry).^2)) + return √(sum((nodes[node_inds[1]].geometry - nodes[node_inds[2]].geometry) .^ 2)) end - function neighbour_inds(ic; nx, ny, leftcut, rightcut, topcut, bottomcut) - ix, iy = ctoxy(ic, nx=nx) + ix, iy = ctoxy(ic; nx=nx) neighbour_x_inds = [] neighbour_y_inds = [] if ix > 1 @@ -624,25 +666,24 @@ function neighbour_inds(ic; nx, ny, leftcut, rightcut, topcut, bottomcut) if iy < ny append!(neighbour_y_inds, iy + 1) end - + neighbour_inds = [] - for x_ind in neighbour_x_inds + for x_ind ∈ neighbour_x_inds append!(neighbour_inds, xytoc(x_ind, iy; nx=nx)) end - for y_ind in neighbour_y_inds + for y_ind ∈ neighbour_y_inds append!(neighbour_inds, xytoc(ix, y_ind; nx=nx)) end return neighbour_inds end - function get_neighbour_inds(ic, gmtry, it) nx = gmtry["dim"]["nx"] ny = gmtry["dim"]["ny"] - ix, iy = ctoxy(ic, nx=nx) + ix, iy = ctoxy(ic; nx=nx) neighbour_inds = [] # println(ix, ", ", iy) - for neighbour in ["left", "right", "top", "bottom"] + for neighbour ∈ ["left", "right", "top", "bottom"] nix = gmtry["data"][neighbour*"ix"][it, iy, ix] + 2 niy = gmtry["data"][neighbour*"iy"][it, iy, ix] + 2 # println(neighbour, ": ", nix, ", ", niy) @@ -653,30 +694,34 @@ function get_neighbour_inds(ic, gmtry, it) return neighbour_inds end - function attach_neightbours(cells, edges, gmtry, it) - for (ic, cell) in enumerate(cells) - for neighbour_ind in get_neighbour_inds(ic, gmtry, it) - for boundary in cell.boundary - for neighbour_boundary in cells[neighbour_ind].boundary - if boundary.index == neighbour_boundary.index && neighbour_ind ∉ boundary.neighbours + for (ic, cell) ∈ enumerate(cells) + for neighbour_ind ∈ get_neighbour_inds(ic, gmtry, it) + for boundary ∈ cell.boundary + for neighbour_boundary ∈ cells[neighbour_ind].boundary + if boundary.index == neighbour_boundary.index && + neighbour_ind ∉ boundary.neighbours append!(boundary.neighbours, neighbour_ind) end end end end end - for (ic, cell) in enumerate(cells) - for edge_ind in [bnd.index for bnd in cell.boundary] - neighbour_edge_inds = [bnd.index for bnd in cell.boundary] - for neighbour_ind in get_neighbour_inds(ic, gmtry, it) - union!(neighbour_edge_inds, [bnd.index for bnd in cells[neighbour_ind].boundary]) + for (ic, cell) ∈ enumerate(cells) + for edge_ind ∈ [bnd.index for bnd ∈ cell.boundary] + neighbour_edge_inds = [bnd.index for bnd ∈ cell.boundary] + for neighbour_ind ∈ get_neighbour_inds(ic, gmtry, it) + union!( + neighbour_edge_inds, + [bnd.index for bnd ∈ cells[neighbour_ind].boundary], + ) end setdiff!(neighbour_edge_inds, edge_ind) - for neighbour_edge_ind in neighbour_edge_inds - for edge_bnd in edges[edge_ind].boundary - for neighbour_edge_bnd in edges[neighbour_edge_ind].boundary - if edge_bnd.index == neighbour_edge_bnd.index && neighbour_edge_ind ∉ edge_bnd.neighbours + for neighbour_edge_ind ∈ neighbour_edge_inds + for edge_bnd ∈ edges[edge_ind].boundary + for neighbour_edge_bnd ∈ edges[neighbour_edge_ind].boundary + if edge_bnd.index == neighbour_edge_bnd.index && + neighbour_edge_ind ∉ edge_bnd.neighbours append!(edge_bnd.neighbours, neighbour_edge_ind) end end @@ -686,14 +731,13 @@ function attach_neightbours(cells, edges, gmtry, it) end end - """ dict2prop!(obj, dict) Copies grid_ggd and space description in dict format to the data structure recursively. """ -function dict2prop!(obj, dict) - for (key, prop) in dict +dict2prop!(obj, dict) = + for (key, prop) ∈ dict if isa(key, Int) if length(obj) < key resize!(obj, key) @@ -707,8 +751,6 @@ function dict2prop!(obj, dict) end end end -end - """ path_to_obj(obj, path) @@ -724,7 +766,7 @@ Note: If integer is -1, the array of field is resized to increase by 1 in length and last element is returned. """ function path_to_obj(obj, path) - for ele in path + for ele ∈ path if typeof(ele) == String obj = getfield(obj, Symbol(ele)) elseif typeof(ele) == Int @@ -740,32 +782,28 @@ function path_to_obj(obj, path) return obj end - isint(x) = typeof(x) == Int - """ val_obj(var, ggd, grid_ggd_index) Given SOLPS variable name (var), return the field to write values on from -ggd object. +ggd object. """ solps_var_to_imas = YAML_load_file("$(@__DIR__)/solps_var_to_imas.yml") -function val_obj(ggd, var, grid_ggd_index) +val_obj(ggd, var, grid_ggd_index) = if var ∉ keys(solps_var_to_imas) return nothing else path = solps_var_to_imas[var] gsi_ind = findlast(isint, path) path_to_prop = path[1:gsi_ind] - prop_to_obj = path[gsi_ind + 1:end] + prop_to_obj = path[gsi_ind+1:end] prop = path_to_obj(ggd, path_to_prop) prop.grid_index = grid_ggd_index prop.grid_subset_index = gsi_ind return path_to_obj(prop, prop_to_obj) end -end - """ find_subset_index() @@ -778,7 +816,7 @@ We'll call them dd_index (the place in the DD) and ggd_index (the type of subset """ function find_subset_index(gsdesc, ggd_index) subsets = length(gsdesc["grid_subset"]) - for dd_index = 1:subsets + for dd_index ∈ 1:subsets if gsdesc["grid_subset"][dd_index]["identifier"]["index"] == ggd_index return dd_index end @@ -786,14 +824,13 @@ function find_subset_index(gsdesc, ggd_index) return 0 # Indicates failure end - """ get_grid_subset_with_index(grid_subset_index) Returns the grid_subset in a grid_ggd with the matching grid_subset_index """ function get_grid_subset_with_index(grid_ggd, grid_subset_index) - for subset in grid_ggd.grid_subset + for subset ∈ grid_ggd.grid_subset if subset.identifier.index == grid_subset_index return subset end @@ -801,24 +838,27 @@ function get_grid_subset_with_index(grid_ggd, grid_subset_index) return 0 # Indicates failure end - -function get_subset_boundary_inds(space::OMAS.edge_profiles__grid_ggd___space, subset::OMAS.edge_profiles__grid_ggd___grid_subset) +function get_subset_boundary_inds( + space::OMAS.edge_profiles__grid_ggd___space, + subset::OMAS.edge_profiles__grid_ggd___grid_subset, +) nD = subset.element[1].object[1].dimension if nD > 0 - nD_objects = space.objects_per_dimension[nD + 1].object - elements = [nD_objects[ele.object[1].index] for ele in subset.element] + nD_objects = space.objects_per_dimension[nD+1].object + elements = [nD_objects[ele.object[1].index] for ele ∈ subset.element] boundary_inds = Int[] - for ele in elements - symdiff!(boundary_inds, [bnd.index for bnd in ele.boundary]) - + for ele ∈ elements + symdiff!(boundary_inds, [bnd.index for bnd ∈ ele.boundary]) end return boundary_inds end return [] end - -function get_subset_boundary(space::OMAS.edge_profiles__grid_ggd___space, subset::OMAS.edge_profiles__grid_ggd___grid_subset) +function get_subset_boundary( + space::OMAS.edge_profiles__grid_ggd___space, + subset::OMAS.edge_profiles__grid_ggd___grid_subset, +) ret_subset = OMAS.edge_profiles__grid_ggd___grid_subset() boundary_inds = get_subset_boundary_inds(space, subset) bnd_dim = subset.element[1].object[1].dimension - 1 @@ -827,33 +867,41 @@ function get_subset_boundary(space::OMAS.edge_profiles__grid_ggd___space, subset return ret_subset.element end - function get_subset_space(space::OMAS.edge_profiles__grid_ggd___space, - elements::Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}) + elements::Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}) nD = elements[1].object[1].dimension nD_objects = space.objects_per_dimension[nD+1].object - return [nD_objects[ele.object[1].index] for ele in elements] + return [nD_objects[ele.object[1].index] for ele ∈ elements] end """ - subset_do(set_operator, itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}}; - space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(), - use_nodes=false) - -Function to perform any set operation (intersect, union, setdiff etc.) on subset.element to -generate a list of elements to go to subset object. If use_nodes is true, the set operation will be -applied on the set of nodes from subset.element, space argument is required for this. -Note: that the arguments are subset.element (not the subset itself). Similarly, the return object is a -list of OMAS.edge_profiles__grid_ggd___grid_subset___element. + subset_do(set_operator, + itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}}; + space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(), + use_nodes=false) + +Function to perform any set operation (intersect, union, setdiff etc.) on +subset.element to generate a list of elements to go to subset object. If use_nodes is +true, the set operation will be applied on the set of nodes from subset.element, space +argument is required for this. +Note: that the arguments are subset.element (not the subset itself). Similarly, the +return object is a list of OMAS.edge_profiles__grid_ggd___grid_subset___element. """ -function subset_do(set_operator, itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}}; - space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(), - use_nodes=false) +function subset_do(set_operator, + itrs::Vararg{Vector{OMAS.edge_profiles__grid_ggd___grid_subset___element}}; + space::OMAS.edge_profiles__grid_ggd___space=OMAS.edge_profiles__grid_ggd___space(), + use_nodes=false) if use_nodes - ele_inds = set_operator([union([obj.nodes for obj in get_subset_space(space, set_elements)]...) for set_elements in itrs]...) + ele_inds = set_operator( + [ + union([obj.nodes for obj ∈ get_subset_space(space, set_elements)]...) for set_elements ∈ itrs + ]..., + ) dim = 0 else - ele_inds = set_operator([[ele.object[1].index for ele in set_elements] for set_elements in itrs]...) + ele_inds = set_operator( + [[ele.object[1].index for ele ∈ set_elements] for set_elements ∈ itrs]..., + ) dim = itrs[1][1].object[1].dimension end ret_subset = OMAS.edge_profiles__grid_ggd___grid_subset() @@ -862,7 +910,6 @@ function subset_do(set_operator, itrs::Vararg{Vector{OMAS.edge_profiles__grid_gg return ret_subset.element end - """ solps2imas(b2gmtry, b2output, gsdesc; load_bb=false) @@ -897,62 +944,61 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) cut_keys = ["leftcut", "rightcut", "bottomcut", "topcut"] cuts_found = cut_keys ⊆ keys(gmtry["data"]) if cuts_found - cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key in cut_keys]) + cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key ∈ cut_keys]) xpoints_nodes = get_xpoint_nodes(gmtry) end - if typeof(gsdesc) == String gsdesc = YAML_load_file(gsdesc) end # Add grid_ggd array equal to number of time steps resize!(ids.edge_profiles.grid_ggd, gmtry["dim"]["time"]) - for it in 1:gmtry["dim"]["time"] - + for it ∈ 1:gmtry["dim"]["time"] grid_ggd = ids.edge_profiles.grid_ggd[it] grid_ggd.time = Float64.(gmtry["data"]["timesa"][it]) dict2prop!(grid_ggd, gsdesc) - for sn in keys(gsdesc["space"]) + for sn ∈ keys(gsdesc["space"]) space = grid_ggd.space[sn] - # Assuming following to be standard for now. We can add this info through YAML as well + # Assuming following to be standard for now. + # We can add this info through YAML as well resize!(space.objects_per_dimension, 4) o0 = space.objects_per_dimension[1] # 0D objects o1 = space.objects_per_dimension[2] # 1D objects o2 = space.objects_per_dimension[3] # 2D objects o3 = space.objects_per_dimension[4] # 3D objects - subset_nodes = get_grid_subset_with_index(grid_ggd, 1) # nodes have index 1 - subset_faces = get_grid_subset_with_index(grid_ggd, 2) # faces (edges with elongation in third dimension) have index 2 - subset_xedges = get_grid_subset_with_index(grid_ggd, 3) # faces (edges with elongation in third dimension) have index 2 - subset_yedges = get_grid_subset_with_index(grid_ggd, 4) # faces (edges with elongation in third dimension) have index 2 - subset_cells = get_grid_subset_with_index(grid_ggd, 5) # cells (cell in 2D grid, volume with elongation) have index 5 + subset_nodes = get_grid_subset_with_index(grid_ggd, 1) + subset_faces = get_grid_subset_with_index(grid_ggd, 2) + subset_xedges = get_grid_subset_with_index(grid_ggd, 3) + subset_yedges = get_grid_subset_with_index(grid_ggd, 4) + subset_cells = get_grid_subset_with_index(grid_ggd, 5) if cuts_found - subset_core = get_grid_subset_with_index(grid_ggd, 22) # core cells have index 22 - subset_sol = get_grid_subset_with_index(grid_ggd, 23) # sol cells have index 23 - subset_odr = get_grid_subset_with_index(grid_ggd, 24) # odr cells have index 24 - subset_idr = get_grid_subset_with_index(grid_ggd, 25) # idr cells have index 25 - subset_xp = get_grid_subset_with_index(grid_ggd, 6) # x_points have index 6 - subset_corecut = get_grid_subset_with_index(grid_ggd, 7) # core_cut have index 7 - subset_pfrcut = get_grid_subset_with_index(grid_ggd, 8) # pfr_cut have index 8 - subset_othroat = get_grid_subset_with_index(grid_ggd, 9) # outer_throat index 9 - subset_ithroat = get_grid_subset_with_index(grid_ggd, 10) # inner_throat index 10 + subset_core = get_grid_subset_with_index(grid_ggd, 22) + subset_sol = get_grid_subset_with_index(grid_ggd, 23) + subset_odr = get_grid_subset_with_index(grid_ggd, 24) + subset_idr = get_grid_subset_with_index(grid_ggd, 25) + subset_xp = get_grid_subset_with_index(grid_ggd, 6) + subset_corecut = get_grid_subset_with_index(grid_ggd, 7) + subset_pfrcut = get_grid_subset_with_index(grid_ggd, 8) + subset_othroat = get_grid_subset_with_index(grid_ggd, 9) + subset_ithroat = get_grid_subset_with_index(grid_ggd, 10) if !isnothing(jxa) - subset_omp = get_grid_subset_with_index(grid_ggd, 11) # outer midplane 11 - subset_ompsep = get_grid_subset_with_index(grid_ggd, 101) # outer midplane node at separatix 101 + subset_omp = get_grid_subset_with_index(grid_ggd, 11) + subset_ompsep = get_grid_subset_with_index(grid_ggd, 101) end if !isnothing(jxi) - subset_imp = get_grid_subset_with_index(grid_ggd, 12) # inner midplane 12 - subset_impsep = get_grid_subset_with_index(grid_ggd, 102) # inner midplane node at separatix 102 + subset_imp = get_grid_subset_with_index(grid_ggd, 12) + subset_impsep = get_grid_subset_with_index(grid_ggd, 102) end - subset_corebnd = get_grid_subset_with_index(grid_ggd, 15) # core inner boundary index 15 - subset_separatix = get_grid_subset_with_index(grid_ggd, 16) # separatix index 16 - subset_otsep = get_grid_subset_with_index(grid_ggd, 103) # outer target separatix meeting point 103 - subset_itsep = get_grid_subset_with_index(grid_ggd, 104) # inner target separatix meeting point 104 + subset_corebnd = get_grid_subset_with_index(grid_ggd, 15) + subset_separatix = get_grid_subset_with_index(grid_ggd, 16) + subset_otsep = get_grid_subset_with_index(grid_ggd, 103) + subset_itsep = get_grid_subset_with_index(grid_ggd, 104) end - - subset_otarget = get_grid_subset_with_index(grid_ggd, 13) # outer_target index 13 - subset_itarget = get_grid_subset_with_index(grid_ggd, 14) # inner_target index 14 + + subset_otarget = get_grid_subset_with_index(grid_ggd, 13) + subset_itarget = get_grid_subset_with_index(grid_ggd, 14) # Resizing objects to hold cell geometry data # Should be fewer than this many points, but this way we won't under-fill @@ -962,19 +1008,19 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) resize!(o3.object, ncell) # Volumes # Initialize geometry for 0D objects(nodes), nodes for 1D objects(edges) - for i = 1:(ncell * 4) + for i ∈ 1:(ncell*4) nodes[i].geometry = [0.0, 0.0] edges[i].nodes = [0, 0] resize!(edges[i].boundary, 2) - for bnd in edges[i].boundary + for bnd ∈ edges[i].boundary bnd.neighbours = Int64[] end end # Initialize nodes and boundaries for cells - for i = 1:(ncell) + for i ∈ 1:(ncell) cells[i].nodes = [0, 0, 0, 0] resize!(cells[i].boundary, 4) - for bnd in cells[i].boundary + for bnd ∈ cells[i].boundary bnd.neighbours = Int64[] end end @@ -982,19 +1028,25 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) j = 1 edge_ind = 1 # Setting up space with nodes, edges and cells - for iy = 1:ny - for ix = 1:nx + for iy ∈ 1:ny + for ix ∈ 1:nx ic::Int = (iy - 1) * nx + ix - # Adding node positions data to grid_ggd[grid_number].space[space_number].objects_per_dimension[0].object[:].geometry - # Adding cell corners data to grid_ggd[grid_number].space[space_number].objects_per_dimension[2].object[:].nodes[1:4] - for icorner = 1:4 - # Have to search to see if the node is already added and then record its index - # If not already listed, then list it under new index and record that - # Note that time index has been fixed to 1 here. Only handling fixed grid geometry - # through the run cases. - i_existing = search_points(ids, crx[1, icorner, iy, ix], cry[1, icorner, iy, ix])[1] + # Adding node positions and cell corners data + for icorner ∈ 1:4 + # Have to search to see if the node is already added and then + # record its index + # If not already listed, then list it under new index and + # record that + # Note that time index has been fixed to 1 here. Only handling + # fixed grid geometry through the run cases. + i_existing = search_points( + ids, + crx[1, icorner, iy, ix], + cry[1, icorner, iy, ix], + )[1] if i_existing == 0 - nodes[j].geometry = [crx[1, icorner, iy, ix], cry[1, icorner, iy, ix]] + nodes[j].geometry = + [crx[1, icorner, iy, ix], cry[1, icorner, iy, ix]] cells[ic].nodes[icorner] = j add_subset_element!(subset_nodes, sn, 0, j) if cuts_found && xpoints_nodes[it][1] == nodes[j].geometry @@ -1005,21 +1057,36 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) cells[ic].nodes[icorner] = i_existing[1] end end - # Adding edges (faced with toroidal elongation) to grid_ggd[grid_number].space[space_number].objects_per_dimension[1].object[:].nodes[1:2] - # Adding same edges as boundary to grid_ggd[grid_number].space[space_number].objects_per_dimension[2].object[:].boundary[1:4] - for (boundary_ind, edge_pair) in chosen_edge_order - edge_nodes = [cells[ic].nodes[icorner] for icorner in edge_pair] + # Adding edges (faces with toroidal elongation) + # Adding same edges as boundary to cells + for (boundary_ind, edge_pair) ∈ chosen_edge_order + edge_nodes = [cells[ic].nodes[icorner] for icorner ∈ edge_pair] existing_edge_ind = search_edges(edges, edge_nodes) if existing_edge_ind == 0 edges[edge_ind].nodes = edge_nodes - for (ii, edge_bnd) in enumerate(edges[edge_ind].boundary) + for (ii, edge_bnd) ∈ enumerate(edges[edge_ind].boundary) edge_bnd.index = edge_nodes[ii] end - edges[edge_ind].measure = distance_between_nodes(nodes, edge_nodes) + edges[edge_ind].measure = + distance_between_nodes(nodes, edge_nodes) cells[ic].boundary[boundary_ind].index = edge_ind add_subset_element!(subset_faces, sn, 1, edge_ind) - add_subset_element!(subset_xedges, sn, 1, edge_ind, is_x_aligned; boundary_ind) - add_subset_element!(subset_yedges, sn, 1, edge_ind, is_y_aligned; boundary_ind) + add_subset_element!( + subset_xedges, + sn, + 1, + edge_ind, + is_x_aligned; + boundary_ind, + ) + add_subset_element!( + subset_yedges, + sn, + 1, + edge_ind, + is_y_aligned; + boundary_ind, + ) edge_ind += 1 else cells[ic].boundary[boundary_ind].index = existing_edge_ind @@ -1027,11 +1094,37 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) end add_subset_element!(subset_cells, sn, 2, ic) if cuts_found - # add_subset_element!(subset, sn, dim, index, ix, iy, in_subset=(x...) -> true; kwargs...) - add_subset_element!(subset_core, sn, 2, ic, in_core; ix, iy, cuts...) + add_subset_element!( + subset_core, + sn, + 2, + ic, + in_core; + ix, + iy, + cuts..., + ) add_subset_element!(subset_sol, sn, 2, ic, in_sol; iy, cuts...) - add_subset_element!(subset_idr, sn, 2, ic, in_idr; ix, iy, cuts...) - add_subset_element!(subset_odr, sn, 2, ic, in_odr; ix, iy, cuts...) + add_subset_element!( + subset_idr, + sn, + 2, + ic, + in_idr; + ix, + iy, + cuts..., + ) + add_subset_element!( + subset_odr, + sn, + 2, + ic, + in_odr; + ix, + iy, + cuts..., + ) end end end @@ -1039,24 +1132,93 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) # Add boundaries attach_neightbours(cells, edges, gmtry, it) # Adding edges to subsets - for iy = 1:ny - for ix = 1:nx - for boundary_ind = 1:4 - edge_ind = cells[xytoc(ix, iy; nx)].boundary[boundary_ind].index - add_subset_element!(subset_corecut, sn, 1, edge_ind, is_core_cut; ix, iy, cells, nx, boundary_ind, cuts...) - # add_subset_element!(subset_pfrcut, sn, 1, edge_ind, is_pfr_cut; ix, iy, cells, nx, boundary_ind, cuts...) + for iy ∈ 1:ny + for ix ∈ 1:nx + for boundary_ind ∈ 1:4 + edge_ind = + cells[xytoc(ix, iy; nx)].boundary[boundary_ind].index + add_subset_element!( + subset_corecut, + sn, + 1, + edge_ind, + is_core_cut; + ix, + iy, + cells, + nx, + boundary_ind, + cuts..., + ) if !isnothing(jxa) - add_subset_element!(subset_omp, sn, 1, edge_ind, is_outer_midplane; ix, iy, jxa, boundary_ind, cuts...) + add_subset_element!( + subset_omp, + sn, + 1, + edge_ind, + is_outer_midplane; + ix, + iy, + jxa, + boundary_ind, + cuts..., + ) end if !isnothing(jxi) - add_subset_element!(subset_imp, sn, 1, edge_ind, is_inner_midplane; ix, iy, jxi, boundary_ind, cuts...) + add_subset_element!( + subset_imp, + sn, + 1, + edge_ind, + is_inner_midplane; + ix, + iy, + jxi, + boundary_ind, + cuts..., + ) end - add_subset_element!(subset_othroat, sn, 1, edge_ind, is_outer_throat; ix, iy, boundary_ind, cuts...) - add_subset_element!(subset_ithroat, sn, 1, edge_ind, is_inner_throat; ix, iy, boundary_ind, cuts...) - add_subset_element!(subset_otarget, sn, 1, edge_ind, is_outer_target; ix, nx, boundary_ind) - add_subset_element!(subset_itarget, sn, 1, edge_ind, is_inner_target; ix, boundary_ind) - # add_subset_element!(subset_corebnd, sn, 1, edge_ind, is_core_boundary; ix, iy, boundary_ind, cuts...) - # add_subset_element!(subset_separatix, sn, 1, edge_ind, is_separatix; iy, boundary_ind, cuts...) + add_subset_element!( + subset_othroat, + sn, + 1, + edge_ind, + is_outer_throat; + ix, + iy, + boundary_ind, + cuts..., + ) + add_subset_element!( + subset_ithroat, + sn, + 1, + edge_ind, + is_inner_throat; + ix, + iy, + boundary_ind, + cuts..., + ) + add_subset_element!( + subset_otarget, + sn, + 1, + edge_ind, + is_outer_target; + ix, + nx, + boundary_ind, + ) + add_subset_element!( + subset_itarget, + sn, + 1, + edge_ind, + is_inner_target; + ix, + boundary_ind, + ) end end end @@ -1064,20 +1226,46 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) sol_boundary_elements = get_subset_boundary(space, subset_sol) idr_boundary_elements = get_subset_boundary(space, subset_idr) odr_boundary_elements = get_subset_boundary(space, subset_odr) - subset_pfrcut.element = subset_do(intersect, idr_boundary_elements, odr_boundary_elements) - subset_corebnd.element = subset_do(setdiff, core_boundary_elements, sol_boundary_elements) + subset_pfrcut.element = + subset_do(intersect, idr_boundary_elements, odr_boundary_elements) + subset_corebnd.element = + subset_do(setdiff, core_boundary_elements, sol_boundary_elements) subset_separatix.element = subset_do(intersect, sol_boundary_elements, - subset_do(union, core_boundary_elements, - odr_boundary_elements, - idr_boundary_elements)) + subset_do(union, core_boundary_elements, + odr_boundary_elements, + idr_boundary_elements)) if !isnothing(jxa) - subset_ompsep.element = subset_do(intersect, subset_separatix.element, subset_omp.element; space, use_nodes=true) + subset_ompsep.element = subset_do( + intersect, + subset_separatix.element, + subset_omp.element; + space, + use_nodes=true, + ) end if !isnothing(jxi) - subset_impsep.element = subset_do(intersect, subset_separatix.element, subset_imp.element; space, use_nodes=true) + subset_impsep.element = subset_do( + intersect, + subset_separatix.element, + subset_imp.element; + space, + use_nodes=true, + ) end - subset_otsep.element = subset_do(intersect, subset_separatix.element, subset_otarget.element; space, use_nodes=true) - subset_itsep.element = subset_do(intersect, subset_separatix.element, subset_itarget.element; space, use_nodes=true) + subset_otsep.element = subset_do( + intersect, + subset_separatix.element, + subset_otarget.element; + space, + use_nodes=true, + ) + subset_itsep.element = subset_do( + intersect, + subset_separatix.element, + subset_itarget.element; + space, + use_nodes=true, + ) end end # End of setting up space end @@ -1086,15 +1274,15 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) b2 = read_b2_output(b2output) # Add grid_ggd array equal to number of time steps resize!(ids.edge_profiles.ggd, b2["dim"]["time"]) - for it in 1:b2["dim"]["time"] + for it ∈ 1:b2["dim"]["time"] ggd = ids.edge_profiles.ggd[it] ggd.time = Float64.(b2["data"]["timesa"][it]) - for (key, data) in b2["data"] + for (key, data) ∈ b2["data"] obj = val_obj(ggd, key, gsdesc["identifier"]["index"]) if !isnothing(obj) resize!(obj, ncell) - for iy = 1:ny - for ix = 1:nx + for iy ∈ 1:ny + for ix ∈ 1:nx ic::Int = (iy - 1) * nx + ix obj[ic] = data[it, iy, ix] end @@ -1111,7 +1299,7 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) resize!(ids.equilibrium.time_slice, gmtry["dim"]["time"]) ids.equilibrium.time = gmtry["data"]["timesa"] end - for it in 1:gmtry["dim"]["time"] + for it ∈ 1:gmtry["dim"]["time"] # Note # Ideally, equilibrium keeps separate grids_ggd object for each time step # But since we have already created them in edge_profiles.grid_ggd, we @@ -1128,13 +1316,14 @@ function solps2imas(b2gmtry, b2output, gsdesc, b2mn=nothing; load_bb=false) b_r = time_slice.ggd[1].b_field_r[1] b_z = time_slice.ggd[1].b_field_z[1] - b_z.grid_index = b_r.grid_index = b_t.grid_index = gsdesc["identifier"]["index"] + b_z.grid_index = + b_r.grid_index = b_t.grid_index = gsdesc["identifier"]["index"] b_z.grid_subset_index = b_r.grid_subset_index = b_t.grid_subset_index = 5 resize!(b_z.values, ncell) resize!(b_r.values, ncell) resize!(b_t.values, ncell) - for iy = 1:ny - for ix = 1:nx + for iy ∈ 1:ny + for ix ∈ 1:nx ic::Int = (iy - 1) * nx + ix b_z.values[ic] = bb[it, 1, iy, ix] b_r.values[ic] = bb[it, 2, iy, ix] diff --git a/test/runtests.jl b/test/runtests.jl index 46278c6..c4e0f12 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -import SOLPS2IMAS +using SOLPS2IMAS: SOLPS2IMAS using Test using YAML: load_file as YAML_load_file @@ -12,18 +12,27 @@ function test_ind_conversion() nx = 92 ny = 38 success = true - for iy = 1:ny - for ix = 1:nx + for iy ∈ 1:ny + for ix ∈ 1:nx cxy = SOLPS2IMAS.ctoxy(SOLPS2IMAS.xytoc(ix, iy; nx=nx); nx=nx) if cxy != (ix, iy) ic = SOLPS2IMAS.xytoc(ix, iy; nx=nx) - println("ic: ", ic, ", (ix, iy): (", ix, ", ", iy, "), converted_ix,iy: ", cxy) + println( + "ic: ", + ic, + ", (ix, iy): (", + ix, + ", ", + iy, + "), converted_ix,iy: ", + cxy, + ) success = false end end end - for ic = 1:nx*ny + for ic ∈ 1:nx*ny cic = SOLPS2IMAS.xytoc(SOLPS2IMAS.ctoxy(ic; nx=nx)...; nx=nx) if cic != ic ix, iy = SOLPS2IMAS.ctoxy(ic; nx=nx) @@ -69,7 +78,7 @@ function test_read_b2_output() @assert(size(contents["data"]["fl3dr"]) == (nt, nybr)) @assert(size(contents["data"]["fc3dr"]) == (nt, nybr)) @assert(size(contents["data"]["fna3da"]) == (nt, ns, nybr)) - @assert(size(contents["data"]["tmhacore"]) == (nt, )) + @assert(size(contents["data"]["tmhacore"]) == (nt,)) @assert(size(contents["data"]["tmhasol"]) == (nt,)) @assert(size(contents["data"]["tmhadiv"]) == (nt,)) return true @@ -88,29 +97,37 @@ function test_solps2imas() it = 3 iy = 4 ix = 5 - @assert(b2t["data"]["ne2d"][3, iy, ix] == dd.edge_profiles.ggd[3].electrons.density[5].values[(iy - 1) * nx + ix]) + @assert( + b2t["data"]["ne2d"][3, iy, ix] == + dd.edge_profiles.ggd[3].electrons.density[5].values[(iy-1)*nx+ix] + ) # Checking if subsets obtained from set operations are same # as using a brute force definition which is too dependent # on the correct ordering of nodes in SOLPS data files. cut_keys = ["leftcut", "rightcut", "bottomcut", "topcut"] gmtry = SOLPS2IMAS.read_b2_output(b2gmtry) ny = gmtry["dim"]["ny"] - cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key in cut_keys]) - subset_pfrcut = SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 8) - subset_corebnd = SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 15) - subset_separatix= SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 16) + cuts = Dict([(Symbol(key), gmtry["data"][key][1]) for key ∈ cut_keys]) + subset_pfrcut = + SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 8) + subset_corebnd = + SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 15) + subset_separatix = + SOLPS2IMAS.get_grid_subset_with_index(dd.edge_profiles.grid_ggd[1], 16) cells = dd.edge_profiles.grid_ggd[1].space[1].objects_per_dimension[3].object - subset_pfrcut_element_list = [ele.object[1].index for ele in subset_pfrcut.element] - subset_corebnd_element_list = [ele.object[1].index for ele in subset_corebnd.element] - subset_separatix_element_list = [ele.object[1].index for ele in subset_separatix.element] + subset_pfrcut_element_list = [ele.object[1].index for ele ∈ subset_pfrcut.element] + subset_corebnd_element_list = [ele.object[1].index for ele ∈ subset_corebnd.element] + subset_separatix_element_list = + [ele.object[1].index for ele ∈ subset_separatix.element] brute_force_pfrcut_list = [] brute_force_corebnd_list = [] brute_force_separatix_list = [] - for iy = 1:ny - for ix = 1:nx - for boundary_ind = 1:4 - edge_ind = cells[SOLPS2IMAS.xytoc(ix, iy; nx)].boundary[boundary_ind].index - if SOLPS2IMAS.is_pfr_cut(;ix, iy, cells, nx, boundary_ind, cuts...) + for iy ∈ 1:ny + for ix ∈ 1:nx + for boundary_ind ∈ 1:4 + edge_ind = + cells[SOLPS2IMAS.xytoc(ix, iy; nx)].boundary[boundary_ind].index + if SOLPS2IMAS.is_pfr_cut(; ix, iy, cells, nx, boundary_ind, cuts...) append!(brute_force_pfrcut_list, edge_ind) elseif SOLPS2IMAS.is_core_boundary(; ix, iy, boundary_ind, cuts...) append!(brute_force_corebnd_list, edge_ind) @@ -126,7 +143,6 @@ function test_solps2imas() return true end - @testset "omasstuff" begin @test SOLPS2IMAS.try_omas() === nothing @test test_ind_conversion()