diff --git a/Examples/parastell_example.py b/Examples/parastell_example.py index 65d678a9..ba32b6ef 100644 --- a/Examples/parastell_example.py +++ b/Examples/parastell_example.py @@ -40,7 +40,7 @@ "shield": {"thickness_matrix": uniform_unit_thickness * 50}, "vacuum_vessel": { "thickness_matrix": uniform_unit_thickness * 10, - "h5m_tag": "vac_vessel", + "mat_tag": "vac_vessel", }, } # Construct in-vessel components diff --git a/Examples/radial_build_distance_example.py b/Examples/radial_build_distance_example.py deleted file mode 100644 index ded2c9f3..00000000 --- a/Examples/radial_build_distance_example.py +++ /dev/null @@ -1,36 +0,0 @@ -from parastell.radial_distance_utils import * - -coils_file = "coils.example" -width = 40.0 -thickness = 50.0 -toroidal_extent = 90.0 -vmec_file = "wout_vmec.nc" -vmec = read_vmec.VMECData(vmec_file) - -magnet_set = magnet_coils.MagnetSet( - coils_file, width, thickness, toroidal_extent -) - -reordered_filaments = get_reordered_filaments(magnet_set) - -build_magnet_surface(reordered_filaments) - -toroidal_angles = np.linspace(0, 90, num=4) -poloidal_angles = np.linspace(0, 360, num=4) -wall_s = 1.08 - -radial_build_dict = {} - -radial_build = ivb.RadialBuild( - toroidal_angles, - poloidal_angles, - wall_s, - radial_build_dict, - split_chamber=False, -) -build = ivb.InVesselBuild(vmec, radial_build, num_ribs=72, num_rib_pts=96) -build.populate_surfaces() -build.calculate_loci() -ribs = build.Surfaces["chamber"].Ribs -radial_distances = measure_radial_distance(ribs) -np.savetxt("radial_distances.csv", radial_distances, delimiter=",") diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py new file mode 100644 index 00000000..22faf13f --- /dev/null +++ b/Examples/radial_distance_example.py @@ -0,0 +1,114 @@ +import numpy as np + +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 +export_dir = "" +# Define plasma equilibrium VMEC file +vmec_file = "wout_vmec.nc" + +# Instantiate ParaStell build +stellarator = ps.Stellarator(vmec_file) + +# Define build parameters for in-vessel components +# Use 13 x 13 uniformly spaced grid for in-vessel build +toroidal_angles = np.linspace(0, 90, num=13) +poloidal_angles = np.linspace(0, 360, num=13) +wall_s = 1.08 +# Define build parameters for magnet coils +coils_file = "coils.example" +width = 40.0 +thickness = 50.0 +toroidal_extent = 90.0 + +# Measure separation between first wall and coils +available_space = rdu.measure_fw_coils_separation( + vmec_file, + toroidal_angles, + poloidal_angles, + wall_s, + coils_file, + width, + thickness, + sample_mod=1, +) +# 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 = enforce_helical_symmetry(available_space) +# Smooth matrix +steps = 25 # Apply Gaussian filter 25 times +sigma = 0.5 # Smooth using half a standard deviation for the Gaussian kernel +available_space = smooth_matrix(available_space, steps, sigma) +# 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 = enforce_helical_symmetry(available_space) +# Modify available space to account for thickness of magnets +available_space = available_space - max(width, thickness) + +# Define a matrix of uniform unit thickness +uniform_unit_thickness = np.ones((len(toroidal_angles), len(poloidal_angles))) +# Define thickness matrices for each in-vessel component of uniform thickness +first_wall_thickness_matrix = uniform_unit_thickness * 5 +back_wall_thickness_matrix = uniform_unit_thickness * 5 +shield_thickness_matrix = uniform_unit_thickness * 35 +vacuum_vessel_thickness_matrix = uniform_unit_thickness * 30 + +# Compute breeder thickness matrix +breeder_thickness_matrix = ( + available_space + - first_wall_thickness_matrix + - back_wall_thickness_matrix + - shield_thickness_matrix + - vacuum_vessel_thickness_matrix +) + +radial_build_dict = { + "first_wall": {"thickness_matrix": first_wall_thickness_matrix}, + "breeder": {"thickness_matrix": breeder_thickness_matrix}, + "back_wall": {"thickness_matrix": back_wall_thickness_matrix}, + "shield": {"thickness_matrix": shield_thickness_matrix}, + "vacuum_vessel": { + "thickness_matrix": vacuum_vessel_thickness_matrix, + "mat_tag": "vac_vessel", + }, +} + +# Construct in-vessel components +stellarator.construct_invessel_build( + toroidal_angles, + poloidal_angles, + wall_s, + radial_build_dict, + # Set num_ribs and num_rib_pts such that four (60/12 - 1) additional + # locations are interpolated between each specified location + num_ribs=61, + num_rib_pts=61, +) +# Export in-vessel component files +stellarator.export_invessel_build() + +# Construct magnets +stellarator.construct_magnets( + coils_file, width, thickness, toroidal_extent, sample_mod=6 +) +# Export magnet files +stellarator.export_magnets() + +# Define source mesh parameters +mesh_size = (11, 81, 61) +toroidal_extent = 90.0 +# Construct source +stellarator.construct_source_mesh(mesh_size, toroidal_extent) +# Export source file +stellarator.export_source_mesh(filename="source_mesh", export_dir=export_dir) + +# Build Cubit model of Parastell Components +stellarator.build_cubit_model(skip_imprint=False, legacy_faceting=True) + +# Export DAGMC neutronics H5M file +stellarator.export_dagmc(filename="dagmc", export_dir=export_dir) 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 2be72f1d..8a256873 100644 --- a/parastell/magnet_coils.py +++ b/parastell/magnet_coils.py @@ -8,7 +8,7 @@ from . import log from . import cubit_io as cubit_io -from .utils import read_yaml_config, filter_kwargs, m2cm +from .utils import read_yaml_config, filter_kwargs, reorder_loop, m2cm export_allowed_kwargs = ["step_filename", "export_mesh", "mesh_filename"] @@ -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_list = np.array([coil.center_of_mass for coil in filtered_coils]) - com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0]) - # 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 = [ - coil - for _, coil in sorted(zip(com_toroidal_angles, filtered_coils)) - ] + 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,36 +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 create_magnet(self): """Creates a single magnet coil CAD solid in CadQuery. @@ -436,6 +409,83 @@ 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) + shifted_coords = np.append(self.coords[1:], [self.coords[1]], axis=0) + midplane_flags = -np.sign(self.coords[:, 2] * shifted_coords[:, 2]) + # Find index of outboard midplane point + outboard_index = np.argmax(midplane_flags * radii) + + return outboard_index + + def reorder_coords(self, index): + """Reorders coil filament coordinate loop about a given index. + + Arguments: + index (int): index about which to reorder coordinate loop. + """ + self.coords = reorder_loop(self.coords, index) + + def orient_coords(self, positive=True): + """Orients coil filament coordinate loop such that they initially + progress positively or negatively. + + Arguments: + positive (bool): progress coordinates in positive direciton + (defaults to True). If negative, coordinates will progress in + negative direction. + """ + if positive == (self.coords[0, 2] > self.coords[1, 2]): + self.coords = np.flip(self.coords, axis=0) + def parse_args(): """Parser for running as a script""" diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index 69a768ed..35c1ff8b 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -1,154 +1,237 @@ import numpy as np +import cubit +import pystell.read_vmec as read_vmec + from . import magnet_coils from . import invessel_build as ivb -import pystell.read_vmec as read_vmec +from . import cubit_io +from .utils import downsample_loop -import cubit - -def calc_z_radius(point): - """ - Calculate the distance from the z axis. +def reorder_filament(coil): + """Reorders the filament coordinates of a MagnetCoil class object such that + they begin near the outboard midplane, and initially progress positively in + the z-direction. Arguments: - point (iterable of [x,y,z] float): point to find distance for + coil (object): MagnetCoil class object. + Returns: - (float): distance to z axis. + reordered_coords (2-D iterable of float): reordered list of Cartesian + coordinates defining a MagnetCoil filament. """ - return (point[0] ** 2 + point[1] ** 2) ** 0.5 + # Start the filament at the outboard midplane + outboard_index = coil.get_ob_mp_index() + if outboard_index != 0: + coil.reorder_coords(outboard_index) + # Ensure points initially progress in positive z-direction + coil.orient_coords() -def get_start_index(filament): - """ - Find the index at which the filament crosses the xy plane on the OB - side of the filament +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 + progress positively in the z-direction. Arguments: - filament (list of list of float): List of points defining the filament + magnet_set (object): MagnetSet class object. Returns: - crossing_index (int): index at which the filament crosses the xy plane - on the OB side of the filament. + magnet_coils (list of object): reordered list of MagnetCoil class + objects. """ - crossing_index = None - max_crossing_radius = 0 - for index, (point, next_point) in enumerate( - zip(filament[0:-1], filament[1:]) - ): - if point[2] / next_point[2] < 0: - crossing_radius = calc_z_radius(point) - if max_crossing_radius < crossing_radius: - crossing_index = index - max_crossing_radius = crossing_radius - return crossing_index - - -def sort_filaments_toroidally(filaments): - """ - Reorder filaments in order of increasing toroidal angle + magnet_set.populate_magnet_coils() + magnet_coils = magnet_set.magnet_coils + + [reorder_filament(coil) for coil in magnet_coils] + + return magnet_coils + + +def build_line(point_1, point_2): + """Builds a line between two points in Coreform Cubit. Arguments: - filaments (list of list of list of float): List of filaments, which are - lists of points defining each filament. + point_1 (1-D iterable of float): Cartesian coordinates of first point. + point_2 (1-D iterable of float): Cartesian coordinates of second point. + Returns: - filaments (list of list of list of float): filaments in order of - increasing toroidal angle. + curve_id (int): index of curve created in Coreform Cubit. """ - com_list = np.zeros((len(filaments), 3)) + point_1 = " ".join(str(val) for val in point_1) + point_2 = " ".join(str(val) for val in point_2) + cubit.cmd(f"create curve location {point_1} location {point_2}") + curve_id = cubit.get_last_id("curve") - for idx, fil in enumerate(filaments): - com_list[idx] = np.average(fil, axis=0) + return curve_id - phi_arr = np.arctan2(com_list[:, 1], com_list[:, 0]) - filaments = np.array([x for _, x in sorted(zip(phi_arr, filaments))]) +def build_magnet_surface(magnet_coils, sample_mod=1): + """Builds a surface in Coreform Cubit spanning a list of coil filaments. - return filaments + Arguments: + magnet_coils (list of object): list of MagnetCoil class objects, + ordered toroidally. Each MagnetCoil object must also have its + filament coordinates ordered poloidally (see reorder_coils + function). + sample_mod (int): sampling modifier for filament points (defaults to + 1). For a user-defined value n, every nth point will be sampled. + """ + cubit_io.init_cubit() + + surface_lines = [ + [ + build_line(coord, next_coord) + for coord, next_coord in zip( + downsample_loop(coil.coords, sample_mod), + downsample_loop(next_coil.coords, sample_mod), + ) + ] + for coil, next_coil in zip(magnet_coils[:-1], magnet_coils[1:]) + ] + surface_lines = np.array(surface_lines) + + surface_sections = np.reshape( + surface_lines, + ( + len(magnet_coils) - 1, + len(downsample_loop(magnet_coils[0].coords, sample_mod)), + ), + ) + [ + [ + cubit.cmd(f"create surface skin curve {line} {next_line}") + for line, next_line in zip(lines[:-1], lines[1:]) + ] + for lines in surface_sections + ] -def reorder_filaments(filaments): - """ - Reorder the filaments so they start near the outboard xy plane crossing, - and begin by increasing z value. + +def fire_ray(point, direction): + """Performs a ray-firing operation in Coreform Cubit from a reference point + and along a reference direction to determine the distance from that point + to the magnet surface. Arguments: - filaments (list of list of list of float): List of filaments, which are - lists of points defining each filament. + point (iterable of float): Cartesian coordinates of reference point. + direction (iterable of float): reference direction in + Cartesian-coordinate system. + Returns: - filaments (list of list of list of float): Reordered list of filaments, - suitable for building the magnet surface. + distance (float): distance between reference point and magnet surface, + along reference direction. """ - for filament_index, filament in enumerate(filaments): - # start the filament at the outboard - max_z_index = get_start_index(filament) - reordered_filament = np.concatenate( - [filament[max_z_index:], filament[1:max_z_index]] - ) - - # make sure z is increasing initially - if reordered_filament[0, 2] > reordered_filament[1, 2]: - reordered_filament = np.flip(reordered_filament, axis=0) + cubit.cmd(f"create vertex {point[0]} {point[1]} {point[2]}") + vertex_id = cubit.get_last_id("vertex") + + cubit.cmd( + f"create curve location at vertex {vertex_id} " + f"location fire ray location at vertex {vertex_id} " + f"direction {direction[0]} {direction[1]} {direction[2]} at " + "surface all maximum hits 1" + ) - # ensure filament is a closed loop - if max_z_index != 0: - reordered_filament = np.concatenate( - [reordered_filament, [reordered_filament[0]]] - ) + curve_id = cubit.get_last_id("curve") + distance = cubit.get_curve_length(curve_id) - filaments[filament_index] = reordered_filament + return distance - filaments = sort_filaments_toroidally(filaments) - return filaments +def measure_surface_coils_separation(surface): + """Determines the distance between a given Surface class object and a + surface spanning a set of MagnetCoil class objects via a ray-firing + operation in Coreform Cubit. + Arguments: + surface (object): Surface class object. -def get_reordered_filaments(magnet_set): + Returns: + distance_matrix (2-D np.array of float): matrix of distances between + points defining surface.Ribs and magnet surface, along + Rib._normals. """ - Convenience function to get the reordered filament data from a magnet + distance_matrix = np.array( + [ + [ + fire_ray(point, direction) + for point, direction in zip(rib.rib_loci, rib._normals()) + ] + for rib in surface.Ribs + ] + ) + + return distance_matrix + + +def measure_fw_coils_separation( + vmec_file, + toroidal_angles, + poloidal_angles, + wall_s, + coils_file, + width, + thickness, + sample_mod=1, + custom_fw_profile=None, +): + """Measures the distance between a given first wall definition and a set of + magnet filaments, at specified angular locations and along the profile + normal at those angular locations, using ray-firing in Coreform Cubit. + + Arguments: + vmec_file (str): path to plasma equilibrium VMEC file. + toroidal_angles (array of float): toroidal angles at which distances + should be computed [deg]. + poloidal_angles (array of float): poloidal angles at which distances + should be computed [deg]. + wall_s (float): closed flux surface label extrapolation at wall. + coils_file (str): path to coil filament data file. + width (float): width of coil cross-section in toroidal direction [cm]. + thickness (float): thickness of coil cross-section in radial direction + [cm]. + sample_mod (int): sampling modifier for filament points (defaults to + 1). For a user-defined value n, every nth point will be sampled. + custom_fw_profile (2-D iterable of float): thickness matrix defining + first wall profile (defaults to None). + + Returns: + radial_distance_matrix (2-D np.array of float): """ - magnet_set.populate_magnet_coils() + if custom_fw_profile is None: + custom_fw_profile = np.zeros( + (len(toroidal_angles), len(poloidal_angles)) + ) - filtered_filaments = np.array( - [coil.coords for coil in magnet_set.magnet_coils] + vmec = read_vmec.VMECData(vmec_file) + radial_build_dict = {"chamber": {"thickness_matrix": custom_fw_profile}} + + radial_build = ivb.RadialBuild( + toroidal_angles, + poloidal_angles, + wall_s, + radial_build_dict, + split_chamber=False, + ) + invessel_build = ivb.InVesselBuild( + vmec, + radial_build, + # Set num_ribs and num_rib_pts to be less than length of corresponding + # array to ensure that only defined angular locations are used + num_ribs=len(toroidal_angles) - 1, + num_rib_pts=len(poloidal_angles) - 1, ) - filaments = reorder_filaments(filtered_filaments) + invessel_build.populate_surfaces() + invessel_build.calculate_loci() + surface = invessel_build.Surfaces["chamber"] - return filaments + magnet_set = magnet_coils.MagnetSet( + coils_file, width, thickness, toroidal_angles[-1] - toroidal_angles[0] + ) + reordered_coils = reorder_coils(magnet_set) + build_magnet_surface(reordered_coils, sample_mod=sample_mod) -def build_magnet_surface(reordered_filaments): - loops = [] - for fil1, fil2 in zip(reordered_filaments[0:-1], reordered_filaments[1:]): - for index, _ in enumerate(fil1): - loc1 = " ".join(str(val) for val in fil1[index, :]) - loc2 = " ".join(str(val) for val in fil2[index, :]) - cubit.cmd(f"create curve location {loc1} location {loc2}") - loops.append(cubit.get_last_id("curve")) + radial_distance_matrix = measure_surface_coils_separation(surface) - loops = np.array(loops) - loops = np.reshape( - loops, (len(reordered_filaments) - 1, len(reordered_filaments[0])) - ) - for loop in loops: - for line in loop[0:-1]: - cubit.cmd(f"create surface skin curve {line} {line + 1}") - - -def measure_radial_distance(ribs): - distances = [] - for rib in ribs: - distance_subset = [] - for point, direction in zip(rib.rib_loci, rib._normals()): - cubit.cmd(f"create vertex {point[0]} {point[1]} {point[2]}") - vertex_id = cubit.get_last_id("vertex") - cubit.cmd( - f"create curve location at vertex {vertex_id} " - f"location fire ray location at vertex {vertex_id} " - f"direction {direction[0]} {direction[1]} {direction[2]} at " - "surface all maximum hits 1" - ) - curve_id = cubit.get_last_id("curve") - distance = cubit.get_curve_length(curve_id) - distance_subset.append(distance) - distances.append(distance_subset) - return np.array(distances) + return radial_distance_matrix diff --git a/parastell/utils.py b/parastell/utils.py index fcfed50a..1a367cd1 100644 --- a/parastell/utils.py +++ b/parastell/utils.py @@ -1,70 +1,104 @@ import yaml -import math import numpy as np +import math +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 several qualities about the input matrix: + - Regular spacing between angles defining matrix + - Matrix represents a full stellarator period + - The first row corresponds to a toroidal angle at the beginning of a + period (i.e., poloidal symmetry is expected) 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, num_columns = matrix.shape - ang_list_exp = [] + # Ensure rows represent closed loops + matrix[:, -1] = matrix[:, 0] - init_ang = ang_list[0] - final_ang = ang_list[-1] - ang_extent = final_ang - init_ang + # Ensure poloidal symmetry at beginning of period + matrix[0] = np.concatenate( + [ + # Ceil and floor ensure middle element of odd sized array is + # included only once + matrix[0, : math.ceil(num_columns / 2)], + np.flip(matrix[0, : math.floor(num_columns / 2)]), + ] + ) - ang_diff_avg = ang_extent / (num_ang - 1) + # Ensure helical symmetry toroidally and poloidally by mirroring the period + # about both matrix axes + flattened_matrix = matrix.flatten() + flattened_length = len(flattened_matrix) - for ang, next_ang in zip(ang_list[:-1], ang_list[1:]): - n_ang = math.ceil((next_ang - ang) / ang_diff_avg) + first_half = flattened_matrix[: math.ceil(flattened_length / 2)] + last_half = np.flip(flattened_matrix[: math.floor(flattened_length / 2)]) + flattened_matrix = np.concatenate([first_half, last_half]) - ang_list_exp = np.append( - ang_list_exp, np.linspace(ang, next_ang, num=n_ang + 1)[:-1] - ) + matrix = flattened_matrix.reshape((num_rows, num_columns)) - ang_list_exp = np.append(ang_list_exp, ang_list[-1]) + return matrix - return ang_list_exp +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. -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) + Arguments: + list (iterable of float): list to be expanded. + num (int): desired number of entries in expanded list. - return all_data + Returns: + list_exp (iterable of float): expanded list. + """ + list_exp = [] + + init_entry = list[0] + final_entry = list[-1] + extent = final_entry - init_entry + + avg_diff = extent / (num - 1) + + for entry, next_entry in zip(list[:-1], list[1:]): + num_new_entries = int(round((next_entry - entry) / avg_diff)) + + # Don't append the last entry in the created linspace to avoid adding + # it twice when the next created linspace is appended + list_exp = np.append( + list_exp, + np.linspace(entry, next_entry, num=num_new_entries + 1)[:-1], + ) + + list_exp = np.append(list_exp, final_entry) + + return list_exp def filter_kwargs( @@ -101,3 +135,70 @@ 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: + (iterable): reordered closed loop. + """ + return np.concatenate([list[index:], list[1 : index + 1]]) + + +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_matrix = matrix + + for step in range(steps): + smoothed_matrix = np.minimum( + previous_matrix, + gaussian_filter( + previous_matrix, + sigma=sigma, + mode="wrap", + ), + ) + previous_matrix = smoothed_matrix + + return smoothed_matrix