Skip to content

Commit

Permalink
Improve readability and source mesh connectivity
Browse files Browse the repository at this point in the history
  • Loading branch information
@connoramoreno
connoramoreno committed Oct 21, 2023
1 parent 66c15a4 commit 69b69e9
Show file tree
Hide file tree
Showing 4 changed files with 63 additions and 41 deletions.
37 changes: 19 additions & 18 deletions ExampleScript.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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': {
Expand All @@ -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': {
Expand All @@ -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'
}
}
}
Expand All @@ -93,6 +93,7 @@
'cross_section': ['circle', 20],
'start': 3,
'stop': None,
'sample': 6,
'name': 'magnet_coils',
'h5m_tag': 'magnets'
}
Expand Down
10 changes: 7 additions & 3 deletions magnet_coils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -268,14 +268,15 @@ def create_magnets(filaments, cross_section, 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:
file (str): path to magnet coil data file.
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.
Expand Down Expand Up @@ -312,7 +313,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
Expand All @@ -339,6 +340,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)
}
Expand All @@ -363,7 +367,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
Expand Down
13 changes: 7 additions & 6 deletions parastell.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -495,16 +495,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(
Expand Down Expand Up @@ -581,6 +579,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)
}
Expand Down
44 changes: 30 additions & 14 deletions source_mesh.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -182,25 +182,38 @@ 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.
Returns:
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

Expand Down Expand Up @@ -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, ],
Expand All @@ -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,
Expand All @@ -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, ],
Expand All @@ -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,
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 69b69e9

Please sign in to comment.