From 4f124c6668867fb187d2759f853a9bc4506ee819 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Wed, 18 Sep 2024 12:09:22 -0500 Subject: [PATCH 01/16] Modify code structure --- Examples/parastell_example.py | 2 +- Examples/radial_build_distance_example.py | 178 +++++++++-- parastell/radial_distance_utils.py | 349 +++++++++++++++------- 3 files changed, 404 insertions(+), 125 deletions(-) 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 index ded2c9f3..8d4e4355 100644 --- a/Examples/radial_build_distance_example.py +++ b/Examples/radial_build_distance_example.py @@ -1,36 +1,176 @@ -from parastell.radial_distance_utils import * +import numpy as np +from scipy.ndimage import gaussian_filter +import parastell.parastell as ps +import parastell.radial_distance_utils as rdu + + +def smooth_matrix(matrix, steps): + """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. Analagous to Gaussian sigma. + + 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=1) + ) + previous_iteration = smoothed_matrix + + return smoothed_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 +toroidal_angles = np.linspace(0, 90, num=21) +poloidal_angles = np.linspace(0, 360, num=23) +wall_s = 1.08 + +# Define build parameters for magnet coils 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 +# 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=2, ) +# Modify available space to account for thickness of magnets +available_space = available_space - np.sqrt(2) * thickness / 2 -reordered_filaments = get_reordered_filaments(magnet_set) +# Ensure poloidal symmetry at toroidal angles 0 and 45 degrees +for index in range( + len(poloidal_angles) - 1, int((len(poloidal_angles) - 1) / 2), -1 +): + available_space[0, index] = np.flip( + available_space[0, len(poloidal_angles) - 1 - index] + ) + available_space[int((len(toroidal_angles) - 1) / 2), index] = np.flip( + available_space[ + int((len(toroidal_angles) - 1) / 2), + len(poloidal_angles) - 1 - index, + ] + ) +# Ensure quasi-symmetry toroidally and poloidally +for index in range( + len(toroidal_angles) - 1, int((len(toroidal_angles) - 1) / 2), -1 +): + available_space[index] = np.flip( + available_space[len(toroidal_angles) - 1 - index] + ) -build_magnet_surface(reordered_filaments) +available_space = smooth_matrix(available_space, 100) -toroidal_angles = np.linspace(0, 90, num=4) -poloidal_angles = np.linspace(0, 360, num=4) -wall_s = 1.08 +# Ensure poloidal symmetry at toroidal angles 0 and 45 degrees +for index in range( + len(poloidal_angles) - 1, int((len(poloidal_angles) - 1) / 2), -1 +): + available_space[0, index] = np.flip( + available_space[0, len(poloidal_angles) - 1 - index] + ) + available_space[int((len(toroidal_angles) - 1) / 2), index] = np.flip( + available_space[ + int((len(toroidal_angles) - 1) / 2), + len(poloidal_angles) - 1 - index, + ] + ) +# Ensure quasi-symmetry toroidally and poloidally +for index in range( + len(toroidal_angles) - 1, int((len(toroidal_angles) - 1) / 2), -1 +): + available_space[index] = np.flip( + available_space[len(toroidal_angles) - 1 - index] + ) + +print(available_space) -radial_build_dict = {} +# 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 * 50 +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 = ivb.RadialBuild( +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, - split_chamber=False, + num_ribs=61, + num_rib_pts=67, ) -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=",") +# Export in-vessel component files +stellarator.export_invessel_build( + export_cad_to_dagmc=False, export_dir=export_dir +) + +# Construct magnets +stellarator.construct_magnets( + coils_file, width, thickness, toroidal_extent, sample_mod=6 +) +# Export magnet files +stellarator.export_magnets( + step_filename="magnets", + # export_mesh=True, + mesh_filename="magnet_mesh", + export_dir=export_dir, +) +""" +# 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/radial_distance_utils.py b/parastell/radial_distance_utils.py index 69a768ed..70f33b1d 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -1,154 +1,293 @@ 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 -import cubit +def calc_radius(point): + """Calculates the cylindrical radius of a reference point. + + Arguments: + point (iterable of float): Cartesian coordinates of reference + point. -def calc_z_radius(point): + Returns: + (float): cylindrical radius. """ - Calculate the distance from the z axis. + return np.linalg.norm(point[0:2]) + + +def get_start_index(coil): + """Finds the index of the outboard midplane coordinate on a coil filament. Arguments: - point (iterable of [x,y,z] float): point to find distance for + coil (object): MagnetCoil class object. + Returns: - (float): distance to z axis. + outboard_index (int): index of the outboard midplane point. """ - return (point[0] ** 2 + point[1] ** 2) ** 0.5 + radii = [calc_radius(point) for point in coil.coords] + # Determine whether adjacent points cross the midplane + midplane_flags = np.less( + coil.coords[:, 2] / np.append(coil.coords[1:, 2], coil.coords[1, 2]), + np.zeros(len(coil.coords)), + ) + # Find index of outboard midplane point + outboard_index = np.argmax(midplane_flags * radii) + return outboard_index -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_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: - filament (list of list of float): List of points defining the filament + coil (object): MagnetCoil class object. Returns: - crossing_index (int): index at which the filament crosses the xy plane - on the OB side of the filament. + reordered_coords (2-D iterable of float): reordered list of Cartesian + coordinates defining a MagnetCoil filament. """ - 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 + # 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]] + ) + # Ensure filament is a closed loop + if outboard_index != 0: + reordered_coords = np.concatenate( + [reordered_coords, [reordered_coords[0]]] + ) + # 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: - filaments (list of list of list of float): List of filaments, which are - lists of points defining each filament. + magnet_coils (list of object): list of MagnetCoil class objects. + Returns: - filaments (list of list of list of float): filaments in order of - increasing toroidal angle. + magnet_coils (list of object): sorted list of MagnetCoil class objects. """ - com_list = np.zeros((len(filaments), 3)) + com_list = np.array([coil.center_of_mass for coil in magnet_coils]) + com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0]) - for idx, fil in enumerate(filaments): - com_list[idx] = np.average(fil, axis=0) + magnet_coils = np.array( + [x for _, x in sorted(zip(com_toroidal_angles, magnet_coils))] + ) - phi_arr = np.arctan2(com_list[:, 1], com_list[:, 0]) + return magnet_coils - filaments = np.array([x for _, x in sorted(zip(phi_arr, filaments))]) - return filaments +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: + magnet_set (object): MagnetSet class object. -def reorder_filaments(filaments): + Returns: + magnet_coils (list of object): reordered list of MagnetCoil class + objects. """ - Reorder the filaments so they start near the outboard xy plane crossing, - and begin by increasing z value. + magnet_set.populate_magnet_coils() + magnet_coils = magnet_set.magnet_coils + + [reorder_filament(coil) for coil in magnet_coils] + magnet_coils = sort_coils_toroidally(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): Reordered list of filaments, - suitable for building the magnet surface. + curve_id (int): index of curve created in Coreform Cubit. """ - 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]] - ) + 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") + + return curve_id + - # make sure z is increasing initially - if reordered_filament[0, 2] > reordered_filament[1, 2]: - reordered_filament = np.flip(reordered_filament, axis=0) +def build_magnet_surface(magnet_coils, sample_mod=1): + """Builds a surface in Coreform Cubit spanning a list of coil filaments. - # ensure filament is a closed loop - if max_z_index != 0: - reordered_filament = np.concatenate( - [reordered_filament, [reordered_filament[0]]] + 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( + np.concatenate( + [coil.coords[:-1:sample_mod], [coil.coords[0]]] + ), + np.concatenate( + [next_coil.coords[:-1:sample_mod], [next_coil.coords[0]]] + ), ) + ] + 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(magnet_coils[0].coords[:-1:sample_mod]) + 1, + ), + ) - filaments[filament_index] = reordered_filament + [ + [ + 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 + ] - filaments = sort_filaments_toroidally(filaments) - return filaments +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: + point (iterable of float): Cartesian coordinates of reference point. + direction (iterable of float): reference direction in + Cartesian-coordinate system. -def get_reordered_filaments(magnet_set): + Returns: + distance (float): distance between reference point and magnet surface, + along reference direction. """ - Convenience function to get the reordered filament data from a magnet + 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) + + return distance + + +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. + + Returns: + distance_matrix (2-D np.array of float): matrix of distances between + points defining surface.Ribs and magnet surface, along + Rib._normals. """ - magnet_set.populate_magnet_coils() + distance_matrix = [ + [ + fire_ray(point, distance) + for point, distance in zip(rib.rib_loci, rib._normals()) + ] + for rib in surface.Ribs + ] + + return np.array(distance_matrix) + + +def measure_fw_coils_separation( + vmec_file, + toroidal_angles, + poloidal_angles, + wall_s, + coils_file, + width, + thickness, + sample_mod=1, +): + """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. - filtered_filaments = np.array( - [coil.coords for coil in magnet_set.magnet_coils] + Returns: + radial_distance_matrix (2-D np.array of float): + """ + vmec = read_vmec.VMECData(vmec_file) + radial_build_dict = {} + + 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 From fc183b01d0f28c8f52ddcf98e8132df18d059465 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 19 Sep 2024 11:27:15 -0500 Subject: [PATCH 02/16] Add utility functions --- Examples/radial_build_distance_example.py | 109 +++++----------------- parastell/radial_distance_utils.py | 79 ++++++++++++++-- 2 files changed, 95 insertions(+), 93 deletions(-) diff --git a/Examples/radial_build_distance_example.py b/Examples/radial_build_distance_example.py index 8d4e4355..fa4364a2 100644 --- a/Examples/radial_build_distance_example.py +++ b/Examples/radial_build_distance_example.py @@ -1,32 +1,9 @@ import numpy as np -from scipy.ndimage import gaussian_filter import parastell.parastell as ps import parastell.radial_distance_utils as rdu -def smooth_matrix(matrix, steps): - """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. Analagous to Gaussian sigma. - - 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=1) - ) - previous_iteration = smoothed_matrix - - return smoothed_matrix - - # Define directory to export all output files to export_dir = "" # Define plasma equilibrium VMEC file @@ -36,10 +13,9 @@ def smooth_matrix(matrix, steps): stellarator = ps.Stellarator(vmec_file) # Define build parameters for in-vessel components -toroidal_angles = np.linspace(0, 90, num=21) -poloidal_angles = np.linspace(0, 360, num=23) +toroidal_angles = np.linspace(0, 90, num=61) +poloidal_angles = np.linspace(0, 360, num=67) wall_s = 1.08 - # Define build parameters for magnet coils coils_file = "coils.example" width = 40.0 @@ -55,64 +31,26 @@ def smooth_matrix(matrix, steps): coils_file, width, thickness, - sample_mod=2, + sample_mod=1, ) +# For matrices defined by angles that are regularly spaced, measurement results +# in matrix elements that are close to, but not exactly, helcially symmetric +available_space = rdu.enforce_helical_symmetry(available_space) +# Smooth matrix +available_space = rdu.smooth_matrix(available_space, 50, 1) +# For matrices defined by angles that are regularly spaced, matrix smoothing +# results in matrix elements that are close to, but not exactly, helcially +# symmetric +available_space = rdu.enforce_helical_symmetry(available_space) # Modify available space to account for thickness of magnets -available_space = available_space - np.sqrt(2) * thickness / 2 - -# Ensure poloidal symmetry at toroidal angles 0 and 45 degrees -for index in range( - len(poloidal_angles) - 1, int((len(poloidal_angles) - 1) / 2), -1 -): - available_space[0, index] = np.flip( - available_space[0, len(poloidal_angles) - 1 - index] - ) - available_space[int((len(toroidal_angles) - 1) / 2), index] = np.flip( - available_space[ - int((len(toroidal_angles) - 1) / 2), - len(poloidal_angles) - 1 - index, - ] - ) -# Ensure quasi-symmetry toroidally and poloidally -for index in range( - len(toroidal_angles) - 1, int((len(toroidal_angles) - 1) / 2), -1 -): - available_space[index] = np.flip( - available_space[len(toroidal_angles) - 1 - index] - ) - -available_space = smooth_matrix(available_space, 100) - -# Ensure poloidal symmetry at toroidal angles 0 and 45 degrees -for index in range( - len(poloidal_angles) - 1, int((len(poloidal_angles) - 1) / 2), -1 -): - available_space[0, index] = np.flip( - available_space[0, len(poloidal_angles) - 1 - index] - ) - available_space[int((len(toroidal_angles) - 1) / 2), index] = np.flip( - available_space[ - int((len(toroidal_angles) - 1) / 2), - len(poloidal_angles) - 1 - index, - ] - ) -# Ensure quasi-symmetry toroidally and poloidally -for index in range( - len(toroidal_angles) - 1, int((len(toroidal_angles) - 1) / 2), -1 -): - available_space[index] = np.flip( - available_space[len(toroidal_angles) - 1 - index] - ) - -print(available_space) +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 * 50 +shield_thickness_matrix = uniform_unit_thickness * 35 vacuum_vessel_thickness_matrix = uniform_unit_thickness * 30 # Compute breeder thickness matrix @@ -134,31 +72,28 @@ def smooth_matrix(matrix, steps): "mat_tag": "vac_vessel", }, } +# radial_build_dict = {"space": {"thickness_matrix": available_space}} + # Construct in-vessel components stellarator.construct_invessel_build( toroidal_angles, poloidal_angles, wall_s, radial_build_dict, - num_ribs=61, - num_rib_pts=67, + # 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, ) # Export in-vessel component files -stellarator.export_invessel_build( - export_cad_to_dagmc=False, export_dir=export_dir -) +stellarator.export_invessel_build() # Construct magnets stellarator.construct_magnets( coils_file, width, thickness, toroidal_extent, sample_mod=6 ) # Export magnet files -stellarator.export_magnets( - step_filename="magnets", - # export_mesh=True, - mesh_filename="magnet_mesh", - export_dir=export_dir, -) +stellarator.export_magnets() """ # Define source mesh parameters mesh_size = (11, 81, 61) diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index 70f33b1d..c1f52858 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -1,6 +1,7 @@ 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 @@ -217,15 +218,17 @@ def measure_surface_coils_separation(surface): points defining surface.Ribs and magnet surface, along Rib._normals. """ - distance_matrix = [ + distance_matrix = np.array( [ - fire_ray(point, distance) - for point, distance in zip(rib.rib_loci, rib._normals()) + [ + fire_ray(point, distance) + for point, distance in zip(rib.rib_loci, rib._normals()) + ] + for rib in surface.Ribs ] - for rib in surface.Ribs - ] + ) - return np.array(distance_matrix) + return distance_matrix def measure_fw_coils_separation( @@ -291,3 +294,67 @@ 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 From 08c85fbe5ea285d3c99f8755caf087c0eaf50338 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 19 Sep 2024 11:53:26 -0500 Subject: [PATCH 03/16] Rename example file --- ...adial_build_distance_example.py => radial_distance_example.py} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename Examples/{radial_build_distance_example.py => radial_distance_example.py} (100%) diff --git a/Examples/radial_build_distance_example.py b/Examples/radial_distance_example.py similarity index 100% rename from Examples/radial_build_distance_example.py rename to Examples/radial_distance_example.py From 6de1b7a0fa881468732b55775a3ab0ab8fa897a1 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 19 Sep 2024 11:55:23 -0500 Subject: [PATCH 04/16] Uncomment source mesh and DAGMC file generation in example file --- Examples/radial_distance_example.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py index fa4364a2..d5b134bc 100644 --- a/Examples/radial_distance_example.py +++ b/Examples/radial_distance_example.py @@ -94,7 +94,7 @@ ) # Export magnet files stellarator.export_magnets() -""" + # Define source mesh parameters mesh_size = (11, 81, 61) toroidal_extent = 90.0 @@ -108,4 +108,3 @@ # Export DAGMC neutronics H5M file stellarator.export_dagmc(filename="dagmc", export_dir=export_dir) -""" From 68c96e486da306513f472fef550e36250e85e1a9 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 19 Sep 2024 11:56:10 -0500 Subject: [PATCH 05/16] Remove line in example file from testing --- Examples/radial_distance_example.py | 1 - 1 file changed, 1 deletion(-) diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py index d5b134bc..1764c73b 100644 --- a/Examples/radial_distance_example.py +++ b/Examples/radial_distance_example.py @@ -72,7 +72,6 @@ "mat_tag": "vac_vessel", }, } -# radial_build_dict = {"space": {"thickness_matrix": available_space}} # Construct in-vessel components stellarator.construct_invessel_build( From 5ff0d224f9c18a7db5c3caed8c78aea39f9cc440 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 19 Sep 2024 12:32:33 -0500 Subject: [PATCH 06/16] Change documentation wording in example and modify variable name as suggested by Edgar --- Examples/radial_distance_example.py | 7 ++++--- parastell/radial_distance_utils.py | 4 ++-- 2 files changed, 6 insertions(+), 5 deletions(-) diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py index 1764c73b..c6f16fe9 100644 --- a/Examples/radial_distance_example.py +++ b/Examples/radial_distance_example.py @@ -33,13 +33,14 @@ thickness, sample_mod=1, ) -# For matrices defined by angles that are regularly spaced, measurement results -# in matrix elements that are close to, but not exactly, helcially symmetric +# 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) # Smooth matrix available_space = rdu.smooth_matrix(available_space, 50, 1) # For matrices defined by angles that are regularly spaced, matrix smoothing -# results in matrix elements that are close to, but not exactly, helcially +# can result in matrix elements that are close to, but not exactly, helcially # symmetric available_space = rdu.enforce_helical_symmetry(available_space) # Modify available space to account for thickness of magnets diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index c1f52858..8ba6606a 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -221,8 +221,8 @@ def measure_surface_coils_separation(surface): distance_matrix = np.array( [ [ - fire_ray(point, distance) - for point, distance in zip(rib.rib_loci, rib._normals()) + fire_ray(point, direction) + for point, direction in zip(rib.rib_loci, rib._normals()) ] for rib in surface.Ribs ] From a119a83b064aff902df71cfbecf7fa2ee5ec922d Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Fri, 20 Sep 2024 09:48:22 -0500 Subject: [PATCH 07/16] Add ability to define custom first wall profile --- parastell/radial_distance_utils.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index 8ba6606a..8d7ad787 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -240,6 +240,7 @@ def measure_fw_coils_separation( 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 @@ -258,12 +259,19 @@ def measure_fw_coils_separation( [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): """ + if custom_fw_profile is None: + custom_fw_profile = np.zeros( + (len(toroidal_angles), len(poloidal_angles)) + ) + vmec = read_vmec.VMECData(vmec_file) - radial_build_dict = {} + radial_build_dict = {"chamber": {"thickness_matrix": custom_fw_profile}} radial_build = ivb.RadialBuild( toroidal_angles, From 93ffa4a0378731729197202dcd5575ea5d13d562 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Mon, 23 Sep 2024 14:38:23 -0500 Subject: [PATCH 08/16] Incorporate review suggestions --- parastell/magnet_coils.py | 18 ++++++--- parastell/radial_distance_utils.py | 60 +++++++++++++----------------- 2 files changed, 37 insertions(+), 41 deletions(-) diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py index 2be72f1d..72a4601e 100644 --- a/parastell/magnet_coils.py +++ b/parastell/magnet_coils.py @@ -193,16 +193,16 @@ def _filter_coils(self): ] # 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]) + 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) # 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 = sorted( + filtered_coils, key=lambda x: x.com_toroidal_angle() + ) def _cut_magnets(self): """Cuts the magnets at the planes defining the toriodal extent. @@ -370,6 +370,12 @@ def in_toroidal_extent(self, lower_bound, upper_bound): 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. diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index 8d7ad787..71a6d2fc 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -8,19 +8,6 @@ from . import cubit_io -def calc_radius(point): - """Calculates the cylindrical radius of a reference point. - - Arguments: - point (iterable of float): Cartesian coordinates of reference - point. - - Returns: - (float): cylindrical radius. - """ - return np.linalg.norm(point[0:2]) - - def get_start_index(coil): """Finds the index of the outboard midplane coordinate on a coil filament. @@ -30,11 +17,12 @@ def get_start_index(coil): Returns: outboard_index (int): index of the outboard midplane point. """ - radii = [calc_radius(point) for point in coil.coords] - # Determine whether adjacent points cross the midplane - midplane_flags = np.less( - coil.coords[:, 2] / np.append(coil.coords[1:, 2], coil.coords[1, 2]), - np.zeros(len(coil.coords)), + # 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) @@ -78,16 +66,9 @@ def sort_coils_toroidally(magnet_coils): magnet_coils (list of object): list of MagnetCoil class objects. Returns: - magnet_coils (list of object): sorted list of MagnetCoil class objects. + (list of object): sorted list of MagnetCoil class objects. """ - com_list = np.array([coil.center_of_mass for coil in magnet_coils]) - com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0]) - - magnet_coils = np.array( - [x for _, x in sorted(zip(com_toroidal_angles, magnet_coils))] - ) - - return magnet_coils + return sorted(magnet_coils, key=lambda x: x.com_toroidal_angle()) def reorder_coils(magnet_set): @@ -129,6 +110,19 @@ 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. @@ -146,23 +140,19 @@ def build_magnet_surface(magnet_coils, sample_mod=1): [ build_line(coord, next_coord) for coord, next_coord in zip( - np.concatenate( - [coil.coords[:-1:sample_mod], [coil.coords[0]]] - ), - np.concatenate( - [next_coil.coords[:-1:sample_mod], [next_coil.coords[0]]] - ), + 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(magnet_coils[0].coords[:-1:sample_mod]) + 1, + len(downsample_loop(magnet_coils[0].coords, sample_mod)), ), ) From 98dc878040855f46f7f452637514911068d8f746 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Tue, 1 Oct 2024 10:29:09 -0500 Subject: [PATCH 09/16] 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 From 218346fa8baa58b9d8269eb7201529c67ff770af Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Tue, 1 Oct 2024 11:33:17 -0500 Subject: [PATCH 10/16] Modify example --- Examples/radial_distance_example.py | 20 +++++++++++--------- parastell/magnet_coils.py | 24 +++++++++++++++++++++++- parastell/radial_distance_utils.py | 11 +++-------- 3 files changed, 37 insertions(+), 18 deletions(-) diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py index 1ebd24ab..ac3ee1fa 100644 --- a/Examples/radial_distance_example.py +++ b/Examples/radial_distance_example.py @@ -14,8 +14,8 @@ stellarator = ps.Stellarator(vmec_file) # Define build parameters for in-vessel components -toroidal_angles = np.linspace(0, 90, num=61) -poloidal_angles = np.linspace(0, 360, num=67) +toroidal_angles = np.linspace(0, 90, num=97) +poloidal_angles = np.linspace(0, 360, num=97) wall_s = 1.08 # Define build parameters for magnet coils coils_file = "coils.example" @@ -39,7 +39,7 @@ # symmetric available_space = enforce_helical_symmetry(available_space) # Smooth matrix -available_space = smooth_matrix(available_space, 50, 1) +available_space = smooth_matrix(available_space, 100, 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 @@ -48,7 +48,9 @@ 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))) +uniform_unit_thickness = np.ones( + (len(toroidal_angles[::8]), len(poloidal_angles[::8])) +) # 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 @@ -57,7 +59,7 @@ # Compute breeder thickness matrix breeder_thickness_matrix = ( - available_space + available_space[::8, ::8] - first_wall_thickness_matrix - back_wall_thickness_matrix - shield_thickness_matrix @@ -77,14 +79,14 @@ # Construct in-vessel components stellarator.construct_invessel_build( - toroidal_angles, - poloidal_angles, + toroidal_angles[::8], + poloidal_angles[::8], wall_s, radial_build_dict, # 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, + num_ribs=61, + num_rib_pts=61, ) # Export in-vessel component files stellarator.export_invessel_build() diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py index 526e5138..32097e4c 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"] @@ -468,6 +468,28 @@ def get_ob_mp_index(self): 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 and self.coords[0, 2] > self.coords[1, 2]) or ( + not positive and 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 8c409838..35c1ff8b 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -5,7 +5,7 @@ from . import magnet_coils from . import invessel_build as ivb from . import cubit_io -from .utils import reorder_loop, downsample_loop +from .utils import downsample_loop def reorder_filament(coil): @@ -22,15 +22,10 @@ def reorder_filament(coil): """ # Start the filament at the outboard midplane outboard_index = coil.get_ob_mp_index() - - # Ensure filament is a closed loop if outboard_index != 0: - reordered_coords = reorder_loop(coil.coords, outboard_index) - + coil.reorder_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 + coil.orient_coords() def reorder_coils(magnet_set): From 11858f5841e3e8358e259c7977acd19a5cd172d7 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 3 Oct 2024 13:35:54 -0500 Subject: [PATCH 11/16] Update parastell/magnet_coils.py Co-authored-by: Paul Wilson --- parastell/magnet_coils.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py index 32097e4c..728220f4 100644 --- a/parastell/magnet_coils.py +++ b/parastell/magnet_coils.py @@ -485,9 +485,7 @@ def orient_coords(self, positive=True): (defaults to True). If negative, coordinates will progress in negative direction. """ - if (positive and self.coords[0, 2] > self.coords[1, 2]) or ( - not positive and self.coords[0, 2] < self.coords[1, 2] - ): + if (positive == (self.coords[0, 2] > self.coords[1, 2]) ): self.coords = np.flip(self.coords, axis=0) From 9427be29c79bb229b57e3c89b60ac766cc13aa91 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 3 Oct 2024 13:47:10 -0500 Subject: [PATCH 12/16] Update parastell/magnet_coils.py Co-authored-by: Paul Wilson --- parastell/magnet_coils.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py index 728220f4..3e418f7f 100644 --- a/parastell/magnet_coils.py +++ b/parastell/magnet_coils.py @@ -459,10 +459,8 @@ def get_ob_mp_index(self): 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]) - ) + shifted_coords = np.append(self.coords[1:, :], self.coords[1, :]) + midplane_flags = -np.sign(self.coords[:, 2] * shifted_coords[:, 2]) # Find index of outboard midplane point outboard_index = np.argmax(midplane_flags * radii) From 052f2264bbe3bf71bf915636369e0b4a503fe814 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 3 Oct 2024 13:57:42 -0500 Subject: [PATCH 13/16] Update parastell/utils.py Co-authored-by: Paul Wilson --- parastell/utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/parastell/utils.py b/parastell/utils.py index 1718fc08..2caf3bc9 100644 --- a/parastell/utils.py +++ b/parastell/utils.py @@ -161,8 +161,7 @@ def reorder_loop(list, index): 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) + reordered_list = np.concatenate([list[index:], list[1:index+1]]) return reordered_list From 263d3e4640e00eba4489d557c3ba4624da7a7fb2 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 3 Oct 2024 14:03:29 -0500 Subject: [PATCH 14/16] Update parastell/utils.py Co-authored-by: Paul Wilson --- parastell/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parastell/utils.py b/parastell/utils.py index 2caf3bc9..ec9c05f4 100644 --- a/parastell/utils.py +++ b/parastell/utils.py @@ -84,7 +84,7 @@ def expand_list(list, num): np.linspace(entry, next_entry, num=num_new_entries + 1)[:-1], ) - list_exp = np.append(list_exp, list[-1]) + list_exp = np.append(list_exp, final_entry) return list_exp From 12a3fde7d4f2c0dbbaa1008137d82d9ce59bb037 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 3 Oct 2024 14:05:59 -0500 Subject: [PATCH 15/16] Update parastell/utils.py Co-authored-by: Paul Wilson --- parastell/utils.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/parastell/utils.py b/parastell/utils.py index ec9c05f4..ec63be6a 100644 --- a/parastell/utils.py +++ b/parastell/utils.py @@ -79,6 +79,8 @@ def expand_list(list, num): 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], From b1713bfa06dcd5920e54fc322aaa6edb35b1e443 Mon Sep 17 00:00:00 2001 From: Connor Moreno Date: Thu, 3 Oct 2024 15:36:50 -0500 Subject: [PATCH 16/16] Incorporate review suggestions --- Examples/radial_distance_example.py | 23 +++++++++--------- parastell/magnet_coils.py | 4 ++-- parastell/utils.py | 36 +++++++++++++++-------------- 3 files changed, 33 insertions(+), 30 deletions(-) diff --git a/Examples/radial_distance_example.py b/Examples/radial_distance_example.py index ac3ee1fa..22faf13f 100644 --- a/Examples/radial_distance_example.py +++ b/Examples/radial_distance_example.py @@ -14,8 +14,9 @@ stellarator = ps.Stellarator(vmec_file) # Define build parameters for in-vessel components -toroidal_angles = np.linspace(0, 90, num=97) -poloidal_angles = np.linspace(0, 360, num=97) +# 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" @@ -39,7 +40,9 @@ # symmetric available_space = enforce_helical_symmetry(available_space) # Smooth matrix -available_space = smooth_matrix(available_space, 100, 1) +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 @@ -48,9 +51,7 @@ available_space = available_space - max(width, thickness) # Define a matrix of uniform unit thickness -uniform_unit_thickness = np.ones( - (len(toroidal_angles[::8]), len(poloidal_angles[::8])) -) +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 @@ -59,7 +60,7 @@ # Compute breeder thickness matrix breeder_thickness_matrix = ( - available_space[::8, ::8] + available_space - first_wall_thickness_matrix - back_wall_thickness_matrix - shield_thickness_matrix @@ -79,12 +80,12 @@ # Construct in-vessel components stellarator.construct_invessel_build( - toroidal_angles[::8], - poloidal_angles[::8], + toroidal_angles, + poloidal_angles, wall_s, radial_build_dict, - # Set num_ribs and num_rib_pts to be less than length of corresponding - # array to ensure that only defined angular locations are used + # 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, ) diff --git a/parastell/magnet_coils.py b/parastell/magnet_coils.py index 3e418f7f..8a256873 100644 --- a/parastell/magnet_coils.py +++ b/parastell/magnet_coils.py @@ -459,7 +459,7 @@ def get_ob_mp_index(self): 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, :]) + 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) @@ -483,7 +483,7 @@ def orient_coords(self, positive=True): (defaults to True). If negative, coordinates will progress in negative direction. """ - if (positive == (self.coords[0, 2] > self.coords[1, 2]) ): + if positive == (self.coords[0, 2] > self.coords[1, 2]): self.coords = np.flip(self.coords, axis=0) diff --git a/parastell/utils.py b/parastell/utils.py index ec63be6a..2e32a9fa 100644 --- a/parastell/utils.py +++ b/parastell/utils.py @@ -34,24 +34,28 @@ def enforce_helical_symmetry(matrix): num_rows = matrix.shape[0] num_columns = matrix.shape[1] - # Ensure poloidal symmetry at beginning and middle of period + # Ensure rows represent closed loops + for idx in range(num_rows): + matrix[idx, -1] = matrix[idx, 0] + + # Ensure poloidal symmetry at beginning 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]) + # Ensure helical symmetry toroidally and poloidally by mirroring the period + # about both matrix axes + flattened_matrix = matrix.flatten() + + for idx, value in enumerate( + flattened_matrix[: int(len(flattened_matrix) / 2)], start=1 + ): + flattened_matrix[-idx] = value + + matrix = flattened_matrix.reshape((num_rows, num_columns)) return matrix @@ -79,8 +83,8 @@ def expand_list(list, num): 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 + # 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], @@ -161,11 +165,9 @@ def reorder_loop(list, index): index (int): list index about which to reorder loop. Returns: - reordered_loop (iterable): reordered closed loop. + (iterable): reordered closed loop. """ - reordered_list = np.concatenate([list[index:], list[1:index+1]]) - - return reordered_list + return np.concatenate([list[index:], list[1 : index + 1]]) def smooth_matrix(matrix, steps, sigma):