From 98dc878040855f46f7f452637514911068d8f746 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Tue, 1 Oct 2024 10:29:09 -0500 Subject: [PATCH] Restructure code --- Examples/radial_distance_example.py | 7 +- parastell/invessel_build.py | 10 +- parastell/magnet_coils.py | 118 +++++++++++-------- parastell/radial_distance_utils.py | 126 +------------------- parastell/utils.py | 172 +++++++++++++++++++++------- 5 files changed, 218 insertions(+), 215 deletions(-) diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py index c6f16fe9..1ebd24ab 100644 --- a/Examples/radial_distance_example.py +++ b/Examples/radial_distance_example.py @@ -2,6 +2,7 @@ import parastell.parastell as ps import parastell.radial_distance_utils as rdu +from parastell.utils import enforce_helical_symmetry, smooth_matrix # Define directory to export all output files to @@ -36,13 +37,13 @@ # For matrices defined by angles that are regularly spaced, measurement can # result in matrix elements that are close to, but not exactly, helcially # symmetric -available_space = rdu.enforce_helical_symmetry(available_space) +available_space = enforce_helical_symmetry(available_space) # Smooth matrix -available_space = rdu.smooth_matrix(available_space, 50, 1) +available_space = smooth_matrix(available_space, 50, 1) # For matrices defined by angles that are regularly spaced, matrix smoothing # can result in matrix elements that are close to, but not exactly, helcially # symmetric -available_space = rdu.enforce_helical_symmetry(available_space) +available_space = enforce_helical_symmetry(available_space) # Modify available space to account for thickness of magnets available_space = available_space - max(width, thickness) diff --git a/parastell/invessel_build.py b/parastell/invessel_build.py index 1d535ece..b85b269b 100644 --- a/parastell/invessel_build.py +++ b/parastell/invessel_build.py @@ -12,7 +12,7 @@ from . import log from .utils import ( normalize, - expand_ang_list, + expand_list, read_yaml_config, filter_kwargs, m2cm, @@ -181,11 +181,11 @@ def populate_surfaces(self): "Populating surface objects for in-vessel components..." ) - self._toroidal_angles_exp = expand_ang_list( - self.radial_build.toroidal_angles, self.num_ribs + self._toroidal_angles_exp = np.deg2rad( + expand_list(self.radial_build.toroidal_angles, self.num_ribs) ) - self._poloidal_angles_exp = expand_ang_list( - self.radial_build.poloidal_angles, self.num_rib_pts + self._poloidal_angles_exp = np.deg2rad( + expand_list(self.radial_build.poloidal_angles, self.num_rib_pts) ) offset_mat = np.zeros( diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py index 72a4601e..526e5138 100644 --- a/parastell/magnet_coils.py +++ b/parastell/magnet_coils.py @@ -191,18 +191,10 @@ def _filter_coils(self): for coil in self.magnet_coils if coil.in_toroidal_extent(lower_bound, upper_bound) ] - - # Compute toroidal angles of filament centers of mass - com_toroidal_angles = np.array( - [coil.com_toroidal_angle() for coil in filtered_coils] - ) - # Ensure angles are positive - com_toroidal_angles = (com_toroidal_angles + 2 * np.pi) % (2 * np.pi) + self.magnet_coils = filtered_coils # Sort coils by center-of-mass toroidal angle and overwrite stored list - self.magnet_coils = sorted( - filtered_coils, key=lambda x: x.com_toroidal_angle() - ) + self.magnet_coils = self.sort_coils_toroidally() def _cut_magnets(self): """Cuts the magnets at the planes defining the toriodal extent. @@ -299,6 +291,17 @@ def export_mesh(self, mesh_filename="magnet_mesh", export_dir=""): filename=mesh_filename, export_dir=export_dir ) + def sort_coils_toroidally(self): + """Reorders list of coils by toroidal angle on range [-pi, pi]. + + Arguments: + magnet_coils (list of object): list of MagnetCoil class objects. + + Returns: + (list of object): sorted list of MagnetCoil class objects. + """ + return sorted(self.magnet_coils, key=lambda x: x.com_toroidal_angle()) + class MagnetCoil(object): """An object representing a single modular stellarator magnet coil. @@ -340,42 +343,6 @@ def coords(self, data): tangents / np.linalg.norm(tangents, axis=1)[:, np.newaxis] ) - def in_toroidal_extent(self, lower_bound, upper_bound): - """Determines if the coil lies within a given toroidal angular extent, - based on filament coordinates. - - Arguments: - lower_bound (float): lower bound of toroidal extent [rad]. - upper_bound (float): upper bound of toroidal extent [rad]. - - Returns: - in_toroidal_extent (bool): flag to indicate whether coil lies - within toroidal bounds. - """ - # Compute toroidal angle of each point in filament - toroidal_angles = np.arctan2(self._coords[:, 1], self._coords[:, 0]) - # Ensure angles are positive - toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi) - # Compute bounds of toroidal extent of filament - min_tor_ang = np.min(toroidal_angles) - max_tor_ang = np.max(toroidal_angles) - - # Determine if filament toroidal extent overlaps with that of model - if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or ( - max_tor_ang >= lower_bound or max_tor_ang <= upper_bound - ): - in_toroidal_extent = True - else: - in_toroidal_extent = False - - return in_toroidal_extent - - def com_toroidal_angle(self): - """Computes the toroidal angle of the coil center of mass, based on - filament coordinates. - """ - return np.arctan2(self.center_of_mass[1], self.center_of_mass[0]) - def create_magnet(self): """Creates a single magnet coil CAD solid in CadQuery. @@ -442,6 +409,65 @@ def create_magnet(self): shell = cq.Shell.makeShell(face_list) self.solid = cq.Solid.makeSolid(shell) + def in_toroidal_extent(self, lower_bound, upper_bound): + """Determines if the coil lies within a given toroidal angular extent, + based on filament coordinates. + + Arguments: + lower_bound (float): lower bound of toroidal extent [rad]. + upper_bound (float): upper bound of toroidal extent [rad]. + + Returns: + in_toroidal_extent (bool): flag to indicate whether coil lies + within toroidal bounds. + """ + # Compute toroidal angle of each point in filament + toroidal_angles = np.arctan2(self._coords[:, 1], self._coords[:, 0]) + # Ensure angles are positive + toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi) + # Compute bounds of toroidal extent of filament + min_tor_ang = np.min(toroidal_angles) + max_tor_ang = np.max(toroidal_angles) + + # Determine if filament toroidal extent overlaps with that of model + if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or ( + max_tor_ang >= lower_bound or max_tor_ang <= upper_bound + ): + in_toroidal_extent = True + else: + in_toroidal_extent = False + + return in_toroidal_extent + + def com_toroidal_angle(self): + """Computes the toroidal angle of the coil center of mass, based on + filament coordinates. + + Returns: + (float): toroidal angle of coil center of mass [rad]. + """ + return np.arctan2(self.center_of_mass[1], self.center_of_mass[0]) + + def get_ob_mp_index(self): + """Finds the index of the outboard midplane coordinate on a coil + filament. + + Returns: + outboard_index (int): index of the outboard midplane point. + """ + # Compute radial distance of coordinates from z-axis + radii = np.linalg.norm(self.coords[:, :2], axis=1) + # Determine whether adjacent points cross the midplane (if so, they will + # have opposite signs) + midplane_flags = -np.sign( + self.coords[:, 2] + * np.append(self.coords[1:, 2], self.coords[1, 2]) + ) + # Find index of outboard midplane point + outboard_index = np.argmax(midplane_flags * radii) + + return outboard_index + def parse_args(): """Parser for running as a script""" diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index 71a6d2fc..8c409838 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -1,33 +1,11 @@ import numpy as np import cubit import pystell.read_vmec as read_vmec -from scipy.ndimage import gaussian_filter from . import magnet_coils from . import invessel_build as ivb from . import cubit_io - - -def get_start_index(coil): - """Finds the index of the outboard midplane coordinate on a coil filament. - - Arguments: - coil (object): MagnetCoil class object. - - Returns: - outboard_index (int): index of the outboard midplane point. - """ - # Compute radial distance of coordinates from z-axis - radii = np.linalg.norm(coil.coords[:, :2], axis=1) - # Determine whether adjacent points cross the midplane (if so, they will - # have opposite signs) - midplane_flags = -np.sign( - coil.coords[:, 2] * np.append(coil.coords[1:, 2], coil.coords[1, 2]) - ) - # Find index of outboard midplane point - outboard_index = np.argmax(midplane_flags * radii) - - return outboard_index +from .utils import reorder_loop, downsample_loop def reorder_filament(coil): @@ -43,34 +21,18 @@ def reorder_filament(coil): coordinates defining a MagnetCoil filament. """ # Start the filament at the outboard midplane - outboard_index = get_start_index(coil) - reordered_coords = np.concatenate( - [coil.coords[outboard_index:], coil.coords[1:outboard_index]] - ) + outboard_index = coil.get_ob_mp_index() + # Ensure filament is a closed loop if outboard_index != 0: - reordered_coords = np.concatenate( - [reordered_coords, [reordered_coords[0]]] - ) + reordered_coords = reorder_loop(coil.coords, outboard_index) + # Ensure points initially progress in positive z-direction if reordered_coords[0, 2] > reordered_coords[1, 2]: reordered_coords = np.flip(reordered_coords, axis=0) coil.coords = reordered_coords -def sort_coils_toroidally(magnet_coils): - """Reorders list of coils by toroidal angle on range [-pi, pi] (coils - ordered in MagnetSet class by toroidal angle on range [0, 2*pi]). - - Arguments: - magnet_coils (list of object): list of MagnetCoil class objects. - - Returns: - (list of object): sorted list of MagnetCoil class objects. - """ - return sorted(magnet_coils, key=lambda x: x.com_toroidal_angle()) - - def reorder_coils(magnet_set): """Reorders a set of magnetic coils toroidally and reorders their filament coordinates such that they begin near the outboard midplane, and initially @@ -87,7 +49,6 @@ def reorder_coils(magnet_set): magnet_coils = magnet_set.magnet_coils [reorder_filament(coil) for coil in magnet_coils] - magnet_coils = sort_coils_toroidally(magnet_coils) return magnet_coils @@ -110,19 +71,6 @@ def build_line(point_1, point_2): return curve_id -def downsample_loop(list, sample_mod): - """Downsamples a list representing a closed loop. - - Arguments: - points (iterable): closed loop. - sample_mod (int): sampling modifier. - - Returns: - (iterable): downsampled closed loop - """ - return np.concatenate([list[:-1:sample_mod], [list[0]]]) - - def build_magnet_surface(magnet_coils, sample_mod=1): """Builds a surface in Coreform Cubit spanning a list of coil filaments. @@ -292,67 +240,3 @@ def measure_fw_coils_separation( radial_distance_matrix = measure_surface_coils_separation(surface) return radial_distance_matrix - - -def smooth_matrix(matrix, steps, sigma): - """Smooths a matrix via Gaussian filtering, without allowing matrix - elements to increase in value. - - Arguments: - matrix (2-D iterable of float): matrix to be smoothed. - steps (int): number of smoothing steps. - sigma (float): standard deviation for Gaussian kernel. - - Returns: - smoothed_matrix (2-D iterable of float): smoothed matrix. - """ - previous_iteration = matrix - - for step in range(steps): - smoothed_matrix = np.minimum( - previous_iteration, - gaussian_filter( - previous_iteration, - sigma=sigma, - mode="wrap", - ), - ) - previous_iteration = smoothed_matrix - - return smoothed_matrix - - -def enforce_helical_symmetry(matrix): - """Ensures that a matrix is helically symmetric according to stellarator - geometry by overwriting certain matrix elements. Assumes regular spacing - between angles defining matrix. - - Arguments: - matrix (2-D iterable of float): matrix to be made helically symmetric. - - Returns: - matrix (2-D iterable of float): helically symmetric matrix. - """ - num_rows = matrix.shape[0] - num_columns = matrix.shape[1] - - # Ensure poloidal symmetry at beginning and middle of period - matrix[0] = np.concatenate( - [ - matrix[0, : int((num_columns - 1) / 2) + 1], - np.flip(matrix[0, : int(num_columns / 2)]), - ] - ) - matrix[int((num_rows - 1) / 2)] = np.concatenate( - [ - matrix[int((num_rows - 1) / 2), : int((num_columns - 1) / 2) + 1], - np.flip(matrix[int((num_rows - 1) / 2), : int(num_columns / 2)]), - ] - ) - - # Ensure helical symmetry toroidally and poloidally - for index in range(num_rows - 1, int((num_rows - 1) / 2), -1): - matrix[num_rows - 1 - index, -1] = matrix[num_rows - 1 - index, 0] - matrix[index] = np.flip(matrix[num_rows - 1 - index]) - - return matrix diff --git a/parastell/utils.py b/parastell/utils.py index fcfed50a..1718fc08 100644 --- a/parastell/utils.py +++ b/parastell/utils.py @@ -1,70 +1,92 @@ import yaml -import math import numpy as np +from scipy.ndimage import gaussian_filter m2cm = 100 m3tocm3 = m2cm * m2cm * m2cm -def normalize(vec_list): - """Normalizes a set of vectors. +def downsample_loop(list, sample_mod): + """Downsamples a list representing a closed loop. Arguments: - vec_list (1 or 2D np array): single 1D vector or array of 1D vectors - to be normalized + list (iterable): closed loop. + sample_mod (int): sampling modifier. + Returns: - vec_list (np array of same shape as input): single 1D normalized vector - or array of normalized 1D vectors + (iterable): downsampled closed loop """ - if len(vec_list.shape) == 1: - return vec_list / np.linalg.norm(vec_list) - elif len(vec_list.shape) == 2: - return vec_list / np.linalg.norm(vec_list, axis=1)[:, np.newaxis] - else: - print('Input "vec_list" must be 1-D or 2-D NumPy array') + return np.append(list[:-1:sample_mod], [list[0]], axis=0) -def expand_ang_list(ang_list, num_ang): - """Expands list of angles by linearly interpolating according to specified - number to include in stellarator build. +def enforce_helical_symmetry(matrix): + """Ensures that a matrix is helically symmetric according to stellarator + geometry by overwriting certain matrix elements. Assumes regular spacing + between angles defining matrix. Arguments: - ang_list (list of double): user-supplied list of toroidal or poloidal - angles (rad). - num_ang (int): number of angles to include in stellarator build. + matrix (2-D iterable of float): matrix to be made helically symmetric. Returns: - ang_list_exp (list of double): interpolated list of angles (rad). + matrix (2-D iterable of float): helically symmetric matrix. """ - ang_list = np.deg2rad(ang_list) + num_rows = matrix.shape[0] + num_columns = matrix.shape[1] + + # Ensure poloidal symmetry at beginning and middle of period + matrix[0] = np.concatenate( + [ + matrix[0, : int((num_columns - 1) / 2) + 1], + np.flip(matrix[0, : int(num_columns / 2)]), + ] + ) + matrix[int((num_rows - 1) / 2)] = np.concatenate( + [ + matrix[int((num_rows - 1) / 2), : int((num_columns - 1) / 2) + 1], + np.flip(matrix[int((num_rows - 1) / 2), : int(num_columns / 2)]), + ] + ) + + # Ensure helical symmetry toroidally and poloidally + for index in range(num_rows - 1, int((num_rows - 1) / 2), -1): + matrix[num_rows - 1 - index, -1] = matrix[num_rows - 1 - index, 0] + matrix[index] = np.flip(matrix[num_rows - 1 - index]) + + return matrix + + +def expand_list(list, num): + """Expands a list of ordered floats to a total number of entries by + linearly interpolating between entries, inserting a proportional number of + new entries between original entries. - ang_list_exp = [] - - init_ang = ang_list[0] - final_ang = ang_list[-1] - ang_extent = final_ang - init_ang - - ang_diff_avg = ang_extent / (num_ang - 1) + Arguments: + list (iterable of float): list to be expanded. + num (int): desired number of entries in expanded list. - for ang, next_ang in zip(ang_list[:-1], ang_list[1:]): - n_ang = math.ceil((next_ang - ang) / ang_diff_avg) + Returns: + list_exp (iterable of float): expanded list. + """ + list_exp = [] - ang_list_exp = np.append( - ang_list_exp, np.linspace(ang, next_ang, num=n_ang + 1)[:-1] - ) + init_entry = list[0] + final_entry = list[-1] + extent = final_entry - init_entry - ang_list_exp = np.append(ang_list_exp, ang_list[-1]) + avg_diff = extent / (num - 1) - return ang_list_exp + for entry, next_entry in zip(list[:-1], list[1:]): + num_new_entries = int(round((next_entry - entry) / avg_diff)) + list_exp = np.append( + list_exp, + np.linspace(entry, next_entry, num=num_new_entries + 1)[:-1], + ) -def read_yaml_config(filename): - """Read YAML file describing ParaStell configuration and extract all data.""" - with open(filename) as yaml_file: - all_data = yaml.safe_load(yaml_file) + list_exp = np.append(list_exp, list[-1]) - return all_data + return list_exp def filter_kwargs( @@ -101,3 +123,73 @@ def filter_kwargs( raise e return {name: dict[name] for name in allowed_keys} + + +def normalize(vec_list): + """Normalizes a set of vectors. + + Arguments: + vec_list (1 or 2D np array): single 1D vector or array of 1D vectors + to be normalized + Returns: + vec_list (np array of same shape as input): single 1D normalized vector + or array of normalized 1D vectors + """ + if len(vec_list.shape) == 1: + return vec_list / np.linalg.norm(vec_list) + elif len(vec_list.shape) == 2: + return vec_list / np.linalg.norm(vec_list, axis=1)[:, np.newaxis] + else: + print('Input "vec_list" must be 1-D or 2-D NumPy array') + + +def read_yaml_config(filename): + """Read YAML file describing ParaStell configuration and extract all data.""" + with open(filename) as yaml_file: + all_data = yaml.safe_load(yaml_file) + + return all_data + + +def reorder_loop(list, index): + """Reorders a list representing a closed loop. + + Arguments: + list (iterable): closed loop. + index (int): list index about which to reorder loop. + + Returns: + reordered_loop (iterable): reordered closed loop. + """ + reordered_list = np.concatenate([list[index:], list[1:index]]) + reordered_list = np.append(reordered_list, [reordered_list[0]], axis=0) + + return reordered_list + + +def smooth_matrix(matrix, steps, sigma): + """Smooths a matrix via Gaussian filtering, without allowing matrix + elements to increase in value. + + Arguments: + matrix (2-D iterable of float): matrix to be smoothed. + steps (int): number of smoothing steps. + sigma (float): standard deviation for Gaussian kernel. + + Returns: + smoothed_matrix (2-D iterable of float): smoothed matrix. + """ + previous_iteration = matrix + + for step in range(steps): + smoothed_matrix = np.minimum( + previous_iteration, + gaussian_filter( + previous_iteration, + sigma=sigma, + mode="wrap", + ), + ) + previous_iteration = smoothed_matrix + + return smoothed_matrix