Skip to content

Commit

Permalink
Restructure code
Browse files Browse the repository at this point in the history
  • Loading branch information
@connoramoreno
connoramoreno committed Oct 1, 2024
1 parent 93ffa4a commit 98dc878
Show file tree
Hide file tree
Showing 5 changed files with 218 additions and 215 deletions.
7 changes: 4 additions & 3 deletions Examples/radial_distance_example.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down
10 changes: 5 additions & 5 deletions parastell/invessel_build.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from . import log
from .utils import (
normalize,
expand_ang_list,
expand_list,
read_yaml_config,
filter_kwargs,
m2cm,
Expand Down Expand Up @@ -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(
Expand Down
118 changes: 72 additions & 46 deletions parastell/magnet_coils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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.
Expand Down Expand Up @@ -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"""
Expand Down
126 changes: 5 additions & 121 deletions parastell/radial_distance_utils.py
View file
Original file line number Diff line number Diff line change
@@ -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):
Expand All @@ -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
Expand All @@ -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

Expand All @@ -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.
Expand Down Expand Up @@ -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
Loading

0 comments on commit 98dc878

Please sign in to comment.