From 74d70c09f9e25e706402537c948c82bd46c7607a Mon Sep 17 00:00:00 2001 From: Edgar Date: Thu, 29 Aug 2024 17:23:28 -0500 Subject: [PATCH 1/5] radial build distance finder and example --- Examples/radial_build_distance_example.py | 31 ++++ parastell/radial_distance_utils.py | 171 ++++++++++++++++++++++ 2 files changed, 202 insertions(+) create mode 100644 Examples/radial_build_distance_example.py create mode 100644 parastell/radial_distance_utils.py diff --git a/Examples/radial_build_distance_example.py b/Examples/radial_build_distance_example.py new file mode 100644 index 00000000..b34d91d6 --- /dev/null +++ b/Examples/radial_build_distance_example.py @@ -0,0 +1,31 @@ +from parastell.radial_distance_utils import * + +coils_file = "coils_wistelld.txt" +circ_cross_section = ["circle", 25] +toroidal_extent = 90.0 +vmec_file = "plasma_wistelld.nc" +vmec = read_vmec.VMECData(vmec_file) + +magnet_set = magnet_coils.MagnetSet( + coils_file, circ_cross_section, toroidal_extent +) + +reordered_filaments = get_reordered_filaments(magnet_set) + +build_magnet_surface(reordered_filaments) + +toroidal_angles = np.linspace(0, 90, 72) +poloidal_angles = np.linspace(0, 360, 96) +wall_s = 1.08 + +radial_build_dict = {} + +radial_build = ivb.RadialBuild( + toroidal_angles, poloidal_angles, wall_s, radial_build_dict +) +build = ivb.InVesselBuild(vmec, radial_build, split_chamber=False) +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/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py new file mode 100644 index 00000000..e28897ac --- /dev/null +++ b/parastell/radial_distance_utils.py @@ -0,0 +1,171 @@ +import numpy as np + +# from matplotlib import pyplot as plt +# from matplotlib.ticker import FormatStrFormatter +from . import magnet_coils +from . import invessel_build as ivb +import pystell.read_vmec as read_vmec + +import cubit + + +# function to get reordered filaments +def calc_z_radius(point): + """ + Calculate the distance from the z axis. + + Arguments: + point (iterable of [x,y,z] float): point to find distance for + Returns: + (float): distance to z axis. + """ + return (point[0] ** 2 + point[1] ** 2) ** 0.5 + + +def get_start_index(filament): + """ + Find the index at which the filament crosses the xy plane on the OB + side of the filament + + Arguments: + filament (list of list of float): List of points defining the filament + + Returns: + max_z_index (int): index at which the filament crosses the xy plane + """ + max_z_index = None + max_z_radius = None + for index, point in enumerate(filament[0:-1]): + next_point = filament[index + 1] + if point[2] / next_point[2] < 0: + z_radius = calc_z_radius(point) + if max_z_radius is None: + max_z_index = index + max_z_radius = z_radius + + elif max_z_radius < z_radius: + max_z_index = index + max_z_radius = z_radius + return max_z_index + + +def sort_filaments_toroidally(filaments): + """ + Reorder filaments in order of increasing toroidal angle + + Arguments: + filaments (list of list of list of float): List of filaments, which are + lists of points defining each filament. + Returns: + filaments (list of list of list of float): filaments in order of + increasing toroidal angle. + """ + com_list = [] + + for fil in filaments: + com = np.average(fil, axis=0) + com_list.append(com) + + com_list = np.array(com_list) + phi_arr = np.arctan2(com_list[:, 1], com_list[:, 0]) + + filaments = np.array([x for _, x in sorted(zip(phi_arr, filaments))]) + + return filaments + + +def reorder_filaments(filaments): + """ + Reorder the filaments so they start near the outboard xy plane crossing, + and begin by increasing z value. + + Arguments: + filaments (list of list of list of float): List of filaments, which are + lists of points defining each filament. + Returns: + filaments (list of list of list of float): Reorderd list of filaments, + suitable for building the magnet surface. + """ + 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[0:max_z_index]] + ) + + # make sure z is increasing initially + if filament[max_z_index, 2] > filament[max_z_index + 1, 2]: + reordered_filament = np.flip(reordered_filament, axis=0) + + # remove duplicate point since the start point might be different + _, idx = np.unique(reordered_filament, return_index=True, axis=0) + reordered_filament = reordered_filament[np.sort(idx)] + + # ensure filament is a closed loop + reordered_filament = np.concatenate( + [reordered_filament, [reordered_filament[0]]] + ) + + filaments[filament_index] = reordered_filament + + filaments = sort_filaments_toroidally(filaments) + + return filaments + + +def get_reordered_filaments(magnet_set): + """ + Convenience function to get the reordered filament data from a magnet + """ + magnet_set._extract_filaments() + magnet_set._set_average_radial_distance() + magnet_set._set_filtered_filaments() + + filtered_filaments = magnet_set.filtered_filaments + filaments = reorder_filaments(filtered_filaments) + + return filaments + + +def build_magnet_surface(reordered_filaments): + for fil_index, _ in enumerate(reordered_filaments[0:-1]): + fil1 = reordered_filaments[fil_index] + fil2 = reordered_filaments[fil_index + 1] + for index, _ in enumerate(fil1): + x1 = fil1[index, 0] + x2 = fil2[index, 0] + y1 = fil1[index, 1] + y2 = fil2[index, 1] + z1 = fil1[index, 2] + z2 = fil2[index, 2] + cubit.cmd( + f"create curve location {x1} {y1} {z1} location {x2} {y2} {z2}" + ) + + lines = np.array(cubit.get_entities("curve")) + lines = np.reshape( + lines, (len(reordered_filaments) - 1, len(reordered_filaments[0])) + ) + for loop in lines: + 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 = max(cubit.get_entities("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 = max(cubit.get_entities("curve")) + distance = cubit.get_curve_length(curve_id) + distance_subset.append(distance) + distances.append(distance_subset) + return np.array(distances) From d7d33029bdceab6986cc2eb0ed6acf324cbe419c Mon Sep 17 00:00:00 2001 From: Edgar-21 <84034227+Edgar-21@users.noreply.github.com> Date: Fri, 30 Aug 2024 08:52:58 -0500 Subject: [PATCH 2/5] spelling --- parastell/radial_distance_utils.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index e28897ac..c7688ee6 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -9,7 +9,6 @@ import cubit -# function to get reordered filaments def calc_z_radius(point): """ Calculate the distance from the z axis. @@ -83,7 +82,7 @@ def reorder_filaments(filaments): filaments (list of list of list of float): List of filaments, which are lists of points defining each filament. Returns: - filaments (list of list of list of float): Reorderd list of filaments, + filaments (list of list of list of float): Reordered list of filaments, suitable for building the magnet surface. """ for filament_index, filament in enumerate(filaments): From cbc9ab9d98e1527f97049c811f29684f112bd1b3 Mon Sep 17 00:00:00 2001 From: Edgar Date: Tue, 3 Sep 2024 14:52:50 -0500 Subject: [PATCH 3/5] changes per pr comments --- parastell/radial_distance_utils.py | 48 +++++++++++------------------- 1 file changed, 17 insertions(+), 31 deletions(-) diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index c7688ee6..f5428b36 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -1,7 +1,4 @@ import numpy as np - -# from matplotlib import pyplot as plt -# from matplotlib.ticker import FormatStrFormatter from . import magnet_coils from . import invessel_build as ivb import pystell.read_vmec as read_vmec @@ -33,16 +30,13 @@ def get_start_index(filament): max_z_index (int): index at which the filament crosses the xy plane """ max_z_index = None - max_z_radius = None - for index, point in enumerate(filament[0:-1]): - next_point = filament[index + 1] + max_z_radius = 0 + for index, (point, next_point) in enumerate( + zip(filament[0:-1], filament[1:]) + ): if point[2] / next_point[2] < 0: z_radius = calc_z_radius(point) - if max_z_radius is None: - max_z_index = index - max_z_radius = z_radius - - elif max_z_radius < z_radius: + if max_z_radius < z_radius: max_z_index = index max_z_radius = z_radius return max_z_index @@ -59,13 +53,11 @@ def sort_filaments_toroidally(filaments): filaments (list of list of list of float): filaments in order of increasing toroidal angle. """ - com_list = [] + com_list = np.zeros((len(filaments), 3)) - for fil in filaments: - com = np.average(fil, axis=0) - com_list.append(com) + for idx, fil in enumerate(filaments): + com_list[idx] = np.average(fil, axis=0) - com_list = np.array(com_list) phi_arr = np.arctan2(com_list[:, 1], com_list[:, 0]) filaments = np.array([x for _, x in sorted(zip(phi_arr, filaments))]) @@ -89,17 +81,13 @@ def reorder_filaments(filaments): # start the filament at the outboard max_z_index = get_start_index(filament) reordered_filament = np.concatenate( - [filament[max_z_index:], filament[0:max_z_index]] + [filament[max_z_index:], filament[1:max_z_index]] ) # make sure z is increasing initially - if filament[max_z_index, 2] > filament[max_z_index + 1, 2]: + if reordered_filament[0, 2] > reordered_filament[1, 2]: reordered_filament = np.flip(reordered_filament, axis=0) - # remove duplicate point since the start point might be different - _, idx = np.unique(reordered_filament, return_index=True, axis=0) - reordered_filament = reordered_filament[np.sort(idx)] - # ensure filament is a closed loop reordered_filament = np.concatenate( [reordered_filament, [reordered_filament[0]]] @@ -127,9 +115,7 @@ def get_reordered_filaments(magnet_set): def build_magnet_surface(reordered_filaments): - for fil_index, _ in enumerate(reordered_filaments[0:-1]): - fil1 = reordered_filaments[fil_index] - fil2 = reordered_filaments[fil_index + 1] + for fil1, fil2 in zip(reordered_filaments[0:-1], reordered_filaments[1:]): for index, _ in enumerate(fil1): x1 = fil1[index, 0] x2 = fil2[index, 0] @@ -141,11 +127,11 @@ def build_magnet_surface(reordered_filaments): f"create curve location {x1} {y1} {z1} location {x2} {y2} {z2}" ) - lines = np.array(cubit.get_entities("curve")) - lines = np.reshape( - lines, (len(reordered_filaments) - 1, len(reordered_filaments[0])) + loops = np.array(cubit.get_entities("curve")) + loops = np.reshape( + loops, (len(reordered_filaments) - 1, len(reordered_filaments[0])) ) - for loop in lines: + for loop in loops: for line in loop[0:-1]: cubit.cmd(f"create surface skin curve {line} {line + 1}") @@ -156,14 +142,14 @@ def measure_radial_distance(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 = max(cubit.get_entities("vertex")) + 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 = max(cubit.get_entities("curve")) + curve_id = cubit.get_last_id("curve") distance = cubit.get_curve_length(curve_id) distance_subset.append(distance) distances.append(distance_subset) From 75703821c5e122d33fe75438d069b6794db9fd4a Mon Sep 17 00:00:00 2001 From: Edgar Date: Wed, 4 Sep 2024 11:05:19 -0500 Subject: [PATCH 4/5] get curve by last id --- parastell/radial_distance_utils.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index f5428b36..8e9368d0 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -115,6 +115,7 @@ def get_reordered_filaments(magnet_set): def build_magnet_surface(reordered_filaments): + loops = [] for fil1, fil2 in zip(reordered_filaments[0:-1], reordered_filaments[1:]): for index, _ in enumerate(fil1): x1 = fil1[index, 0] @@ -126,8 +127,9 @@ def build_magnet_surface(reordered_filaments): cubit.cmd( f"create curve location {x1} {y1} {z1} location {x2} {y2} {z2}" ) + loops.append(cubit.get_last_id("curve")) - loops = np.array(cubit.get_entities("curve")) + loops = np.array(loops) loops = np.reshape( loops, (len(reordered_filaments) - 1, len(reordered_filaments[0])) ) From 32f7f373c8a754029d9a68124de370c9a6e87137 Mon Sep 17 00:00:00 2001 From: Edgar Date: Wed, 4 Sep 2024 13:36:48 -0500 Subject: [PATCH 5/5] changes per pr comments --- parastell/radial_distance_utils.py | 29 ++++++++++++----------------- 1 file changed, 12 insertions(+), 17 deletions(-) diff --git a/parastell/radial_distance_utils.py b/parastell/radial_distance_utils.py index 8e9368d0..b85e8785 100644 --- a/parastell/radial_distance_utils.py +++ b/parastell/radial_distance_utils.py @@ -27,19 +27,20 @@ def get_start_index(filament): filament (list of list of float): List of points defining the filament Returns: - max_z_index (int): index at which the filament crosses the xy plane + crossing_index (int): index at which the filament crosses the xy plane + on the OB side of the filament. """ - max_z_index = None - max_z_radius = 0 + 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: - z_radius = calc_z_radius(point) - if max_z_radius < z_radius: - max_z_index = index - max_z_radius = z_radius - return max_z_index + 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): @@ -118,15 +119,9 @@ def build_magnet_surface(reordered_filaments): loops = [] for fil1, fil2 in zip(reordered_filaments[0:-1], reordered_filaments[1:]): for index, _ in enumerate(fil1): - x1 = fil1[index, 0] - x2 = fil2[index, 0] - y1 = fil1[index, 1] - y2 = fil2[index, 1] - z1 = fil1[index, 2] - z2 = fil2[index, 2] - cubit.cmd( - f"create curve location {x1} {y1} {z1} location {x2} {y2} {z2}" - ) + 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")) loops = np.array(loops)