From 3da99ddad837c2deda5f060ecad1fe1497a44e11 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Sat, 21 Oct 2023 14:28:27 -0500 Subject: [PATCH] Improve readability and source mesh connectivity --- ExampleScript.py | 37 +++++++++++++++++++------------------ magnet_coils.py | 10 +++++++--- parastell.py | 13 +++++++------ source_mesh.py | 44 ++++++++++++++++++++++++++++++-------------- 4 files changed, 63 insertions(+), 41 deletions(-) diff --git a/ExampleScript.py b/ExampleScript.py index f89089e3..7b0bdd8c 100644 --- a/ExampleScript.py +++ b/ExampleScript.py @@ -19,13 +19,13 @@ [5, 5, 5, 5, 5, 5, 5, 5, 5] ] }, - 'blanket': { + 'breeder': { 'thickness_matrix': [ - [160, 160, 40, 5, 5, 5, 40, 160, 160], - [40, 40, 20, 20, 20, 20, 30, 40, 40 ], - [45, 45, 40, 10, 10, 10, 40, 45, 45 ], - [40, 40, 30, 20, 20, 20, 20, 40, 40 ], - [160, 160, 40, 5, 5, 5, 40, 160, 160] + [100, 100, 30, 10, 10, 10, 30, 100, 100], + [30, 30, 10, 5, 5, 5, 20, 30, 30], + [25, 25, 25, 5, 5, 5, 25, 25, 25], + [30, 30, 20, 5, 5, 5, 10, 30, 30], + [100, 100, 30, 10, 10, 10, 30, 100, 100] ] }, 'back_wall': { @@ -46,13 +46,13 @@ [25, 25, 25, 25, 25, 25, 25, 25, 25] ] }, - 'coolant_manifolds': { + 'manifolds': { 'thickness_matrix': [ - [35, 35, 15, 5, 5, 5, 15, 35, 35], - [15, 15, 5, 5, 5, 5, 5, 15, 15], - [10, 10, 5, 5, 5, 5, 5, 10, 10], - [15, 15, 5, 5, 5, 5, 5, 15, 15], - [35, 35, 15, 5, 5, 5, 15, 35, 35] + [50, 50, 15, 5, 5, 5, 15, 50, 50], + [20, 20, 5, 5, 5, 5, 15, 20, 20], + [15, 15, 15, 5, 5, 5, 15, 15, 15], + [20, 20, 15, 5, 5, 5, 5, 20, 20], + [50, 50, 15, 5, 5, 5, 15, 50, 50] ] }, 'gap': { @@ -69,13 +69,13 @@ # materials with "vacuum" in the name as void material 'vacuum_vessel': { 'thickness_matrix': [ - [20, 20, 20, 20, 20, 20, 20, 20, 20], - [20, 20, 20, 20, 20, 20, 20, 20, 20], - [20, 20, 20, 20, 20, 20, 20, 20, 20], - [20, 20, 20, 20, 20, 20, 20, 20, 20], - [20, 20, 20, 20, 20, 20, 20, 20, 20] + [15, 15, 15, 15, 15, 15, 15, 15, 15], + [15, 15, 15, 15, 15, 15, 15, 15, 15], + [15, 15, 15, 15, 15, 15, 15, 15, 15], + [15, 15, 15, 15, 15, 15, 15, 15, 15], + [15, 15, 15, 15, 15, 15, 15, 15, 15] ], - 'h5m_tag': 'vv' + 'h5m_tag': 'vac_vessel' } } } @@ -93,6 +93,7 @@ 'cross_section': ['circle', 20], 'start': 3, 'stop': None, + 'sample': 6, 'name': 'magnet_coils', 'h5m_tag': 'magnets', 'meshing': True diff --git a/magnet_coils.py b/magnet_coils.py index f2f74a02..9a412420 100644 --- a/magnet_coils.py +++ b/magnet_coils.py @@ -276,7 +276,7 @@ def create_magnets(filaments, cross_section, meshing, logger): return vol_ids -def extract_filaments(file, start, stop): +def extract_filaments(file, start, stop, sample): """Extracts filament data from magnet coil data file. Arguments: @@ -284,6 +284,7 @@ def extract_filaments(file, start, stop): start (int): index for line in data file where coil data begins. stop (int): index for line in data file where coil data ends (defaults to None). + sample (int): sampling modifier for filament points. Returns: filaments (list of list of list of float): list filament coordinates. @@ -320,7 +321,7 @@ def extract_filaments(file, start, stop): # list if s != 0: # Only store every five points - if i % 5 == 0: + if i % sample == 0: # Append coordinates to list coords.append([x, y, z]) # Otherwise, store filament coordinates but do not append final @@ -347,6 +348,9 @@ def magnet_coils(magnets, logger = None): 'cross_section': coil cross-section definition (list), 'start': starting line index for data in file (int), 'stop': stopping line index for data in file (int), + 'sample': sampling modifier for filament points (int). For a + user-supplied value of n, sample every n points in list of + points in each filament, 'name': name to use for STEP export (str), 'h5m_tag': material tag to use in H5M neutronics model (str) } @@ -371,7 +375,7 @@ def magnet_coils(magnets, logger = None): # Extract filament data filaments = extract_filaments( - magnets['file'], magnets['start'], magnets['stop'] + magnets['file'], magnets['start'], magnets['stop'], magnets['sample'] ) # Generate magnet coil solids diff --git a/parastell.py b/parastell.py index 66d24ff0..cbce7e54 100644 --- a/parastell.py +++ b/parastell.py @@ -514,16 +514,14 @@ def expand_ang(ang_list, num_ang): # Compute total angular extent of supplied list ang_ext = ang_list[-1] - ang_list[0] + + # Compute average distance between angles to include in stellarator build + ang_diff_avg = ang_ext/num_ang # Loop over supplied angles for ang, next_ang in zip(ang_list[:-1], ang_list[1:]): - # Compute difference between current and next angles - ang_diff = next_ang - ang - # Compute ratio of angular difference and total angular extent - ratio = ang_diff/ang_ext - # Compute number of angles to interpolate - n_ang = round(ratio*num_ang) + n_ang = int(np.ceil((next_ang - ang)/ang_diff_avg)) # Interpolate angles and append to storage list ang_list_exp = np.append( @@ -600,6 +598,9 @@ def parastell( 'cross_section': coil cross-section definition (list), 'start': starting line index for data in file (int), 'stop': stopping line index for data in file (int), + 'sample': sampling modifier for filament points (int). For a + user-supplied value of n, sample every n points in list of + points in each filament, 'name': name to use for STEP export (str), 'h5m_tag': material tag to use in H5M neutronics model (str) 'meshing': setting for tetrahedral mesh generation (bool) diff --git a/source_mesh.py b/source_mesh.py index c85224c5..9e5ef817 100644 --- a/source_mesh.py +++ b/source_mesh.py @@ -182,13 +182,14 @@ def create_tets_from_wedge( ) -def get_vertex_id(vert_ids, num_s, num_theta): +def get_vertex_id(vert_ids, num_phi, num_s, num_theta): """Computes vertex index in row-major order as stored by MOAB from three-dimensional n x 3 matrix indices. Arguments: vert_ids (list of int): list of vertex indices expressed in n x 3 matrix order. The list format is: [toroidal angle index, flux surface index, poloidal angle index] + num_phi (int): number of toroidal angles defining mesh. num_s (int): number of closed magnetic flux surfaces defining mesh. num_theta (int): number of poloidal angles defining mesh. @@ -196,11 +197,23 @@ def get_vertex_id(vert_ids, num_s, num_theta): id (int): vertex index in row-major order as stored by MOAB """ # Compute index offset from magnetic axis - ma_offset = vert_ids[0] * (1 + ((num_s - 1) * num_theta)) + if vert_ids[0] == num_phi - 1: + # Ensure logical connectivity at final toroidal angle + ma_offset = 0 + else: + # Otherwise, compute index of magnetic axis + ma_offset = vert_ids[0] * (1 + ((num_s - 1) * (num_theta - 1))) # Compute index offset from closed flux surface - s_offset = vert_ids[1] * num_theta + s_offset = vert_ids[1] * (num_theta - 1) + # Compute index offset from poloidal angle + if vert_ids[2] == num_theta: + # Ensure logical connectivity at final poloidal angle + theta_offset = 1 + else: + # Otherwise, compute index of poloidal angle + theta_offset = vert_ids[2] # Compute vertex index - id = ma_offset + s_offset + vert_ids[2] + id = ma_offset + s_offset + theta_offset return id @@ -236,10 +249,10 @@ def create_mesh( # Create tetrahedra, looping through vertices for i in range(num_phi - 1): # Create tetrahedra for wedges at center of plasma - for k in range(num_theta - 1): + for k in range(1, num_theta): # Define six wedge vertex indices in matrix format as well as index # offset accounting for magnetic axis vertex in the form - # [toroidal index, flux surface index, poloidal index, offset] + # [toroidal index, flux surface index, poloidal index] wedge_id_data = [ [i, 0, 0, ], [i, 0, k, ], @@ -252,7 +265,7 @@ def create_mesh( ids_wedge = [] # Define vertex indices in row-major order for wedges for vertex in wedge_id_data: - ids_wedge += [get_vertex_id(vertex, num_s, num_theta)] + ids_wedge += [get_vertex_id(vertex, num_phi, num_s, num_theta)] # Create tetrahedra from wedge create_tets_from_wedge( mbc, tag_handle, tet_set, mbc_verts, verts, verts_s, ids_wedge, @@ -262,11 +275,11 @@ def create_mesh( # Create tetrahedra for hexahedra beyond center of plasma for j in range(num_s - 2): # Create tetrahedra for current hexahedron - for k in range(num_theta - 1): + for k in range(1, num_theta): # Define eight hexahedron vertex indices in matrix format as # well as index offset accounting for magnetic axis vertex in # the form - # [toroidal index, flux surface index, poloidal index, offset] + # [toroidal index, flux surface index, poloidal index] hex_id_data = [ [i, j, k ], [i, j + 1, k, ], @@ -281,7 +294,9 @@ def create_mesh( ids_hex = [] # Define vertex indices in row-major order for hexahedra for vertex in hex_id_data: - ids_hex += [get_vertex_id(vertex, num_s, num_theta)] + ids_hex += [ + get_vertex_id(vertex, num_phi, num_s, num_theta) + ] # Create tetrahedra from hexahedron create_tets_from_hex( mbc, tag_handle, tet_set, mbc_verts, verts, verts_s, @@ -309,8 +324,8 @@ def create_vertices(vmec, mbc, num_s, num_theta, num_phi): vertex. """ # Compute total number of vertices - num_verts = num_phi * (1 + ((num_s - 1) * num_theta)) - # Initialize n x 3 array of verticex coordinates + num_verts = (num_phi - 1) * (1 + ((num_s - 1) * (num_theta - 1))) + # Initialize n x 3 array of vertex coordinates verts = np.zeros((num_verts, 3)) # Initialize array of vertex closed flux surface indices verts_s = np.zeros(num_verts) @@ -331,7 +346,7 @@ def create_vertices(vmec, mbc, num_s, num_theta, num_phi): # Initialize vertex index i_vert = 0 # Compute vertices in Cartesian space - for i_phi in range(num_phi): + for i_phi in range(num_phi - 1): # Compute toroidal angle for vertex phi = i_phi * delta_phi # Determine vertex at magnetic axis, convertic to cm @@ -346,7 +361,8 @@ def create_vertices(vmec, mbc, num_s, num_theta, num_phi): # Compute closed flux surface index for vertex s = i_s * delta_s - for i_theta in range(num_theta): + # Exclude final poloidal angle as the final and initial are the same + for i_theta in range(num_theta - 1): # Compute poloidal angle for vertex theta = i_theta * delta_theta # Determine vertex at specified toroidal angle, closed flux