diff --git a/Project.toml b/Project.toml index a3e1890..1b5c828 100644 --- a/Project.toml +++ b/Project.toml @@ -8,8 +8,10 @@ ArgParse = "c7e460c6-2fb9-53a9-8c5b-16f535851c63" Contour = "d38c429a-6771-53c6-b99e-75d170b6e991" EFIT = "cda752c5-6b03-55a3-9e33-132a441b0c17" GGDUtils = "b7b5e640-9b39-4803-84eb-376048795def" +GeoInterface = "cf35fbd7-0cd7-5166-be24-54bfbe79505f" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6" +LibGEOS = "a90b1aa1-3769-5649-ba7e-abc5a9d163eb" OMAS = "91cfaa06-6526-4804-8666-b540b3feef2f" PhysicalConstants = "5ad8b20f-a522-5ce9-bfc9-ddf1d5bda6ab" PlotUtils = "995b91a9-d308-5afd-9ec6-746e21dbc043" diff --git a/data/.gitignore b/data/.gitignore deleted file mode 100644 index f6bc585..0000000 --- a/data/.gitignore +++ /dev/null @@ -1,2 +0,0 @@ -*.mesh_ext.json - diff --git a/src/SD4SOLPS.jl b/src/SD4SOLPS.jl index 0f87fe4..ab8a12d 100644 --- a/src/SD4SOLPS.jl +++ b/src/SD4SOLPS.jl @@ -8,6 +8,7 @@ using Interpolations: Interpolations export find_files_in_allowed_folders export geqdsk_to_imas! +export preparation include("$(@__DIR__)/supersize_profile.jl") include("$(@__DIR__)/repair_eq.jl") @@ -168,6 +169,9 @@ function geqdsk_to_imas!(eqdsk_file, dd; time_index=1) resize!(limiter.unit, 1) limiter.unit[1].outline.r = g.rlim limiter.unit[1].outline.z = g.zlim + + # Repair + add_rho_to_equilibrium!(dd) return end diff --git a/src/supersize_profile.jl b/src/supersize_profile.jl index cb75739..d8b3bfc 100644 --- a/src/supersize_profile.jl +++ b/src/supersize_profile.jl @@ -8,8 +8,10 @@ using Interpolations: Interpolations using GGDUtils: GGDUtils, get_grid_subset_with_index, add_subset_element!, get_subset_boundary, project_prop_on_subset!, get_subset_centers -using PolygonOps: PolygonOps +using PolygonOps: PolygonOps, inpolygon using JSON: JSON +using LibGEOS: LibGEOS, readgeom, intersection +using GeoInterface: GeoInterface export extrapolate_core export fill_in_extrapolated_core_profile @@ -671,6 +673,70 @@ function modify_mesh_ext_near_x!( end end +""" + itentify_cutoff_cells() + +Finds cells in the extended mesh that are cut off by the wall. Some cells may +intersect the wall and be divided, and some may be outside the wall completely. +Returns two flags: One is "okay" for cells where the cell area matches the area +of the intersection between the cell and the tokamak interior (to within +some small tolerance), and "border" for cells where the area of intersection is +less than this but more than 0. +dd: the data dictionary +r: the R coordinates of the mesh cells to check +z: the Z coordinates of the mesh cells to check +pfr_transition: the index of the mesh row that should not be connected because + it spans the transition from the PFR to the common flux region of the SOL. +""" +function identify_cutoff_cells( + dd::OMAS.dd, + r::Matrix{Float64}, + z::Matrix{Float64}, + pfr_transition::Int64, +) + npol, nlvl = size(r) + okay = zeros(Bool, (npol - 1, nlvl - 1)) + border = zeros(Bool, (npol - 1, nlvl - 1)) + limiter = dd.wall.description_2d[1].limiter + wall_r = limiter.unit[1].outline.r + wall_z = limiter.unit[1].outline.z + poly = [[wall_r[i], wall_z[i]] for i ∈ 1:length(wall_r)] + poly = append!(poly, [[wall_r[1], wall_z[1]]]) + wall_poly = LibGEOS.Polygon([poly]) + # inpoly = [[inpolygon([r[i, j], z[i, j]], poly) for j ∈ 1:nlvl] for i ∈ 1:npol] + # wl = Matrix(reduce(hcat, inpoly)') .!= 0 + tolerance = 1e-7 # m^2 + + for i ∈ 2:npol + if (i - 1) != pfr_transition + for j ∈ 2:nlvl + cell_r = [r[i, j], r[i, j-1], r[i-1, j-1], r[i-1, j], r[i, j]] + cell_z = [z[i, j], z[i, j-1], z[i-1, j-1], z[i-1, j], z[i, j]] + cell = LibGEOS.Polygon([[[cell_r[k], cell_z[k]] for k in 1:length(cell_r)]]) + # within = [wl[i, j], wl[i-1, j], wl[i-1, j-1], wl[i, j-1]] + # okay[i-1, j-1] = all(within) + # border[i-1, j-1] = (minimum(within) == 0) & (maximum(within) == 1) + cell_area = LibGEOS.area(cell) + inter = LibGEOS.intersection(wall_poly, cell) + inter_area = LibGEOS.area(inter) + okay[i-1, j-1] = abs(inter_area - cell_area) < tolerance + border[i-1, j-1] = 0 < inter_area < (cell_area - tolerance) + # println("i=",i,",j=",j,"okay=",okay[i-1,j-1],",ca=",cell_area,",ia=",inter_area) + + # inter = GeoInterface.coordinates(inter)[1] + + end + end + end + # println("area of wall polygon:", LibGEOS.area(wall_poly)) + + # fz = open("limiter.dat", "w") + # for i ∈ 1:length(wall_r) + # println(fz, wall_r[i], " ", wall_z[i]) + # end + return okay, border +end + #!format off """ record_regular_mesh!() @@ -780,7 +846,7 @@ function record_regular_mesh!( add_subset_element!(all_nodes_sub, space_idx, node_dim, node_idx) # Edges - if (i > 1) & (i != cut) # i-1 to i in the npol direction + if (i > 1) & (i != (cut+1)) # i-1 to i in the npol direction edge_idx1 = edge_start1 + iii * e1_per_i + jj edges[edge_idx1].nodes = [node_idx, node_idx - n_per_i] resize!(edges[edge_idx1].boundary, 2) @@ -804,7 +870,7 @@ function record_regular_mesh!( end # Cells - if (i > 1) & (i != cut) & (j > 1) + if (i > 1) & (i != (cut+1)) & (j > 1) cell_idx = cell_start + iii * c_per_i + jjj cells[cell_idx].nodes = [ node_idx, @@ -825,6 +891,184 @@ function record_regular_mesh!( end end +function trim_ext_mesh_with_wall!( + dd::OMAS.dd; + grid_ggd_idx::Int64=1, + space_idx::Int64=1, +) + limiter = dd.wall.description_2d[1].limiter + wall_r = limiter.unit[1].outline.r + wall_z = limiter.unit[1].outline.z + poly = [[wall_r[i], wall_z[i]] for i ∈ 1:length(wall_r)] + poly = append!(poly, [[wall_r[1], wall_z[1]]]) + wall_poly = LibGEOS.Polygon([poly]) + tolerance = 1e-7 # m^2 + + grid_ggd = dd.edge_profiles.grid_ggd[grid_ggd_idx] + space = grid_ggd.space[space_idx] + + o0 = space.objects_per_dimension[1] # Nodes + o1 = space.objects_per_dimension[2] # Edges + o2 = space.objects_per_dimension[3] # Cells (2D projection) + + ext_nodes_sub = get_grid_subset_with_index(grid_ggd, -201) + ext_edges_sub = get_grid_subset_with_index(grid_ggd, -202) + ext_xedges_sub = get_grid_subset_with_index(grid_ggd, -203) + ext_yedges_sub = get_grid_subset_with_index(grid_ggd, -204) + ext_cells_sub = get_grid_subset_with_index(grid_ggd, -205) + + all_nodes_sub = get_grid_subset_with_index(grid_ggd, 1) + all_edges_sub = get_grid_subset_with_index(grid_ggd, 2) + all_xedges_sub = get_grid_subset_with_index(grid_ggd, 3) + all_yedges_sub = get_grid_subset_with_index(grid_ggd, 4) + all_cells_sub = get_grid_subset_with_index(grid_ggd, 5) + all_yedge_indices = [element.object[1].index for element in all_yedges_sub.element] + + deleted_nodes = [] + deleted_edges = [] + deleted_cells = [] + ext_cells = length(ext_cells_sub.element) + for j ∈ 1:ext_cells + i = ext_cells - (j-1) + cell_idx = ext_cells_sub.element[i].object[1].index + nodes = o2.object[cell_idx].nodes + c = [[o0.object[n].geometry[1], o0.object[n].geometry[2]] for n in nodes] + c = [c; [c[1]]] + cell = LibGEOS.Polygon([c]) + cell_area = LibGEOS.area(cell) + inter = LibGEOS.intersection(wall_poly, cell) + inter_area = LibGEOS.area(inter) + if abs(inter_area - cell_area) > tolerance + # print("cell ", i, " has problem: inter=", inter_area, ",cell=",cell_area) + # Problem! This cell needs modification. + edges = o2.object[ext_cells_sub.element[i].object[1].index].boundary + edges = [edge.index for edge in edges] + yedges = [index for index in edges if index in all_yedge_indices] + if inter_area == 0 + # println(", inter area is 0") + # Easy: just delete this guy AND all his nodes + deleted_edges = [deleted_edges; edges] + deleted_nodes = [deleted_nodes; nodes] + + deleteat!(ext_cells_sub.element, i) + # deleteat!(o2.object, cell_idx) + deleted_cells = [deleted_cells; i] + else + # println("cell area is not 0") + # Hard: cell is partly in bounds + # Find which nodes are out of bounds + # Find which edges are out of bounds or intersecting + # Delete all nodes and edges that aren't entirely in bounds + # Delete the cell + # Determine whether new nodes were already added by a previous cell + # The last npol + 1 cells need to be considered based on the order + # in which cells were added. + # Add new nodes that will be unique + # Determine whether new edges were already added by a previous cell + # Add new edges that will be unique + # Create new cell + # Reference the new boundary and nodes + # p = LibGEOS.Point() + + # -------- okay, maybe that would work, but wouldn't it be easier + # to find nodes that are out of bounds and move them along their + # edges until they are on the wall? Then some of the edges would + # be fixed and there would be no invalid nodes. A complicated wall + # shape might require more nodes and edges to be added. + inter = GeoInterface.coordinates(inter)[1] + for yedge in yedges + print(yedge, ": ") + yedge_nodes = o1.object[yedge].nodes + print(yedge_nodes, ": ") + edge_node_coords = [LibGEOS.Point(o0.object[yedge_nodes[1]].geometry[1], o0.object[yedge_nodes[1]].geometry[2]), LibGEOS.Point(o0.object[yedge_nodes[2]].geometry[1], o0.object[yedge_nodes[2]].geometry[2])] + println(edge_node_coords) + edge_line = LibGEOS.LineString([edge_node_coords]) + intersection_point = LibGEOS.intersection(wall_poly, edge_line) + println(intersection_point) + end + end + + end + end + println("# deleted edges = ", length(deleted_edges)) + sub_edges = [ele.object[1].index for ele in ext_edges_sub.element] + sub_xedges = [ele.object[1].index for ele in ext_xedges_sub.element] + sub_yedges = [ele.object[1].index for ele in ext_yedges_sub.element] + asub_edges = [ele.object[1].index for ele in all_edges_sub.element] + asub_xedges = [ele.object[1].index for ele in all_xedges_sub.element] + asub_yedges = [ele.object[1].index for ele in all_yedges_sub.element] + unique!(sort!(deleted_nodes, rev=true)) + unique!(sort!(deleted_edges, rev=true)) + nae = length(asub_edges) + for j in 1:nae + i = nae - (j-1) + if i <= length(sub_edges) + if sub_edges[i] in deleted_edges + deleteat!(ext_edges_sub.element, i) + end + end + if i <= length(sub_xedges) + if sub_xedges[i] in deleted_edges + deleteat!(ext_xedges_sub.element, i) + end + end + if i <= length(sub_yedges) + if sub_yedges[i] in deleted_edges + deleteat!(ext_yedges_sub.element, i) + end + end + if asub_edges[i] in deleted_edges + deleteat!(all_edges_sub.element, i) + end + if i <= length(asub_xedges) + if asub_xedges[i] in deleted_edges + deleteat!(all_xedges_sub.element, i) + end + end + if i <= length(asub_yedges) + if asub_yedges[i] in deleted_edges + deleteat!(all_yedges_sub.element, i) + end + end + end + sub_nodes = [ele.object[1].index for ele in ext_nodes_sub.element] + asub_nodes = [ele.object[1].index for ele in all_nodes_sub.element] + an = length(asub_nodes) + for j in 1:an + i = an - (j-1) + if asub_nodes[i] in deleted_nodes + deleteat!(all_nodes_sub.element, i) + end + if i <= length(sub_nodes) + if sub_nodes[i] in deleted_nodes + deleteat!(ext_nodes_sub.element, i) + end + end + end + asub_cells = [ele.object[1].index for ele in all_cells_sub.element] + nac = length(asub_cells) + for j in 1:nac + i = nac - (j - 1) + if asub_cells[i] in deleted_cells + deleteat!(all_cells_sub.element, i) + end + end + # println("area of wall polygon:", LibGEOS.area(wall_poly)) + # println("length(o0.object)=",length(o0.object), ", maximum(deleted_nodes)=",maximum(deleted_nodes)) + # for node in deleted_nodes + # deleteat!(o0.object, node) + # end + for edge in deleted_edges + deleteat!(o1.object, edge) + end + + fz = open("limiter.dat", "w") + for i ∈ 1:length(wall_r) + println(fz, wall_r[i], " ", wall_z[i]) + end + return +end + """ convert_filename(filename::String) @@ -979,7 +1223,19 @@ function mesh_extension_sol!( # The PFR flag needs to be respected; pfr nodes don't connect to non-pfr nodes pfr = rzpi.(mesh_r[:, 1], mesh_z[:, 1]) .< 1 pfr_transition = argmax(abs.(diff(pfr))) + # okay, border = identify_cutoff_cells(dd,mesh_r,mesh_z,pfr_transition) + # fr = open("okay.dat", "w") + # fz = open("border.dat", "w") + # for i ∈ 1:(npol-1) + # for j ∈ 1:(nlvl-1) + # print(fr, okay[i, j], " ") + # print(fz, border[i, j], " ") + # end + # println(fr, "") + # println(fz, "") + # end record_regular_mesh!(grid_ggd, space, mesh_r, mesh_z, pfr_transition) + trim_ext_mesh_with_wall!(dd; grid_ggd_idx=grid_ggd_idx, space_idx=space_idx) return mesh_r, mesh_z, pfr_transition end