From 94308e6a08b9782bf9cb071bde83a0ae835a0bd4 Mon Sep 17 00:00:00 2001 From: Anchal Gupta Date: Sun, 17 Mar 2024 20:32:32 -0700 Subject: [PATCH] Tests for tri mesh and using resize! 4 nodes edges Correct chosen_tri_edge_order has been identified from the sample EIRENE files. Added some verififcation tests for the EIRENE triangular mesh reading. Still need to check some mismatch in common faces between EIRENE and SOLPS in the PFR edge region. Also now only as many ndoes and edges are being initialized as required in the space. This gets rid of the extra initialization we have been doing so far. --- src/SOLPS2IMAS.jl | 35 +++++++++++++++----------- test/runtests.jl | 64 +++++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 85 insertions(+), 14 deletions(-) diff --git a/src/SOLPS2IMAS.jl b/src/SOLPS2IMAS.jl index bb135a0..c618855 100644 --- a/src/SOLPS2IMAS.jl +++ b/src/SOLPS2IMAS.jl @@ -91,9 +91,9 @@ chosen_edge_order = [(1, (1, 2)), (4, (3, 1))] # Following convention is used to index the edges of a triangular cell -chosen_tri_edge_order = [(1, (2, 3)), - (2, (3, 1)), - (3, (1, 2))] +chosen_tri_edge_order = [(1, (1, 2)), + (2, (2, 3)), + (3, (1, 3))] """ solps2imas( @@ -206,19 +206,21 @@ function solps2imas( # Resizing objects to hold cell geometry data # Should be fewer than this many points, but this way we won't under-fill - nodes = resize!(o1.object, ncell * 4) # Nodes (1D) - edges = resize!(o2.object, ncell * 4) # Edges (2D) + # nodes = resize!(o1.object, ncell * 4) # Nodes (1D) + # edges = resize!(o2.object, ncell * 4) # Edges (2D) + nodes = o1.object # Nodes (1D) + edges = o2.object # Edges (2D) cells = resize!(o3.object, ncell) # Cells (3D) # Initialize geometry for 1D objects(nodes), nodes for 2D objects(edges) - for i ∈ 1:(ncell*4) - nodes[i].geometry = [0.0, 0.0] - edges[i].nodes = [0, 0] - resize!(edges[i].boundary, 2) - for bnd ∈ edges[i].boundary - bnd.neighbours = Int64[] - end - end + # for i ∈ 1:(ncell*4) + # nodes[i].geometry = [0.0, 0.0] + # edges[i].nodes = [0, 0] + # resize!(edges[i].boundary, 2) + # for bnd ∈ edges[i].boundary + # bnd.neighbours = Int64[] + # end + # end # Initialize nodes and boundaries for cells for i ∈ 1:(ncell) cells[i].nodes = [0, 0, 0, 0] @@ -248,6 +250,7 @@ function solps2imas( cry[1, icorner, iy, ix], )[1] if i_existing == 0 + resize!(nodes, j) nodes[j].geometry = [crx[1, icorner, iy, ix], cry[1, icorner, iy, ix]] cells[ic].nodes[icorner] = j @@ -266,7 +269,9 @@ function solps2imas( edge_nodes = [cells[ic].nodes[icorner] for icorner ∈ edge_pair] existing_edge_ind = search_edges(edges, edge_nodes) if existing_edge_ind == 0 + resize!(edges, edge_ind) edges[edge_ind].nodes = edge_nodes + resize!(edges[edge_ind].boundary, 2) for (ii, edge_bnd) ∈ enumerate(edges[edge_ind].boundary) edge_bnd.index = edge_nodes[ii] end @@ -513,7 +518,6 @@ function solps2imas( subset_cells = get_grid_subset_with_index(grid_ggd, 5) subset_b25nodes = grid_ggd.grid_subset[cur_no_subsets+1] subset_b25faces = grid_ggd.grid_subset[cur_no_subsets+2] - subset_b25cells = grid_ggd.grid_subset[cur_no_subsets+3] subset_trinodes = grid_ggd.grid_subset[cur_no_subsets+4] subset_trifaces = grid_ggd.grid_subset[cur_no_subsets+5] subset_tricells = grid_ggd.grid_subset[cur_no_subsets+6] @@ -588,6 +592,9 @@ function solps2imas( resize!(edges, length(edges) + 1) this_edge_ind = length(edges) edges[this_edge_ind].nodes = tri_edge_nodes + for (ii, edge_bnd) ∈ enumerate(edges[this_edge_ind].boundary) + edge_bnd.index = tri_edge_nodes[ii] + end edges[this_edge_ind].measure = distance_between_nodes(nodes, tri_edge_nodes) add_subset_element!(subset_faces, 1, 2, this_edge_ind) diff --git a/test/runtests.jl b/test/runtests.jl index 180d030..8f41840 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -3,6 +3,7 @@ using Test using YAML: load_file as YAML_load_file using ArgParse: ArgParse import OMAS as IMASDD +import SOLPS2IMAS: get_grid_subset_with_index allowed_rtol = 1e-4 @@ -22,6 +23,9 @@ function parse_commandline() ["--parser"], Dict(:help => "Stress test file parsing (other than b2 output files)", :action => :store_true), + ["--fort"], + Dict(:help => "Test triangular mesh generation from fort files", + :action => :store_true), ) args = ArgParse.parse_args(s) if !any(values(args)) # If no flags are set, run all tests @@ -201,3 +205,63 @@ if args["solps2imas"] @test Set(brute_force_separatix_list) == Set(subset_separatix_element_list) end end + +if args["fort"] + @testset "Test triangular mesh generation from fort files" begin + fort = ( + "$(@__DIR__)/../samples/fort.33", + "$(@__DIR__)/../samples/fort.34", + "$(@__DIR__)/../samples/fort.35") + b2gmtry = "$(@__DIR__)/../samples/b2fgmtry" + b2output = "$(@__DIR__)/../samples/b2time.nc" + b2mn = "$(@__DIR__)/../samples/b2mn.dat" + ids = SOLPS2IMAS.solps2imas(b2gmtry, b2output; b2mn=b2mn, fort=fort) + grid_ggd = ids.edge_profiles.grid_ggd[1] + space = grid_ggd.space[1] + + subset_nodes = get_grid_subset_with_index(grid_ggd, 1) + subset_faces = get_grid_subset_with_index(grid_ggd, 2) + subset_cells = get_grid_subset_with_index(grid_ggd, 5) + subset_b25nodes = get_grid_subset_with_index(grid_ggd, -1) + subset_b25faces = get_grid_subset_with_index(grid_ggd, -2) + subset_b25cells = get_grid_subset_with_index(grid_ggd, -5) + subset_trinodes = get_grid_subset_with_index(grid_ggd, -101) + subset_trifaces = get_grid_subset_with_index(grid_ggd, -102) + subset_tricells = get_grid_subset_with_index(grid_ggd, -105) + subset_comnodes = get_grid_subset_with_index(grid_ggd, -201) + subset_comfaces = get_grid_subset_with_index(grid_ggd, -202) + subset_comcells = get_grid_subset_with_index(grid_ggd, -205) + + nodes = grid_ggd.space[1].objects_per_dimension[1].object + edges = grid_ggd.space[1].objects_per_dimension[2].object + cells = grid_ggd.space[1].objects_per_dimension[3].object + + gmtry = SOLPS2IMAS.read_b2_output(b2gmtry) + nx = gmtry["dim"]["nx"] + ny = gmtry["dim"]["ny"] + extra_nodes = 2 * (nx + ny) + + @test length(subset_nodes.element) == length(nodes) + @test length(subset_faces.element) == length(edges) + @test length(subset_cells.element) == length(cells) + + @test length(subset_b25nodes.element) == + length(subset_comnodes.element) + extra_nodes + + @test length(subset_trinodes.element) + extra_nodes == length(nodes) + + for ele ∈ subset_tricells.element + tricell_ind = ele.object[1].index + for bnd_ind ∈ 1:3 + if length(cells[tricell_ind].boundary[bnd_ind].neighbours) > 0 + common_nodes = intersect( + cells[tricell_ind].nodes, + cells[cells[tricell_ind].boundary[bnd_ind].neighbours[1]].nodes, + ) + vert_no = Tuple(indexin(common_nodes, cells[tricell_ind].nodes)) + @test SOLPS2IMAS.chosen_tri_edge_order[bnd_ind][2] == vert_no + end + end + end + end +end