Merge pull request #158 from svalinn/rdu-modifications

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

Examples/radial_distance_example.py

@@ -0,0 +1,114 @@
+import numpy as np
+
+import parastell.parastell as ps
+import parastell.radial_distance_utils as rdu
+from parastell.utils import enforce_helical_symmetry, smooth_matrix
+
+
+# Define directory to export all output files to
+export_dir = ""
+# Define plasma equilibrium VMEC file
+vmec_file = "wout_vmec.nc"
+
+# Instantiate ParaStell build
+stellarator = ps.Stellarator(vmec_file)
+
+# Define build parameters for in-vessel components
+# Use 13 x 13 uniformly spaced grid for in-vessel build
+toroidal_angles = np.linspace(0, 90, num=13)
+poloidal_angles = np.linspace(0, 360, num=13)
+wall_s = 1.08
+# Define build parameters for magnet coils
+coils_file = "coils.example"
+width = 40.0
+thickness = 50.0
+toroidal_extent = 90.0
+
+# Measure separation between first wall and coils
+available_space = rdu.measure_fw_coils_separation(
+    vmec_file,
+    toroidal_angles,
+    poloidal_angles,
+    wall_s,
+    coils_file,
+    width,
+    thickness,
+    sample_mod=1,
+)
+# For matrices defined by angles that are regularly spaced, measurement can
+# result in matrix elements that are close to, but not exactly, helcially
+# symmetric
+available_space = enforce_helical_symmetry(available_space)
+# Smooth matrix
+steps = 25  # Apply Gaussian filter 25 times
+sigma = 0.5  # Smooth using half a standard deviation for the Gaussian kernel
+available_space = smooth_matrix(available_space, steps, sigma)
+# For matrices defined by angles that are regularly spaced, matrix smoothing
+# can result in matrix elements that are close to, but not exactly, helcially
+# symmetric
+available_space = enforce_helical_symmetry(available_space)
+# Modify available space to account for thickness of magnets
+available_space = available_space - max(width, thickness)
+
+# Define a matrix of uniform unit thickness
+uniform_unit_thickness = np.ones((len(toroidal_angles), len(poloidal_angles)))
+# Define thickness matrices for each in-vessel component of uniform thickness
+first_wall_thickness_matrix = uniform_unit_thickness * 5
+back_wall_thickness_matrix = uniform_unit_thickness * 5
+shield_thickness_matrix = uniform_unit_thickness * 35
+vacuum_vessel_thickness_matrix = uniform_unit_thickness * 30
+
+# Compute breeder thickness matrix
+breeder_thickness_matrix = (
+    available_space
+    - first_wall_thickness_matrix
+    - back_wall_thickness_matrix
+    - shield_thickness_matrix
+    - vacuum_vessel_thickness_matrix
+)
+
+radial_build_dict = {
+    "first_wall": {"thickness_matrix": first_wall_thickness_matrix},
+    "breeder": {"thickness_matrix": breeder_thickness_matrix},
+    "back_wall": {"thickness_matrix": back_wall_thickness_matrix},
+    "shield": {"thickness_matrix": shield_thickness_matrix},
+    "vacuum_vessel": {
+        "thickness_matrix": vacuum_vessel_thickness_matrix,
+        "mat_tag": "vac_vessel",
+    },
+}
+
+# Construct in-vessel components
+stellarator.construct_invessel_build(
+    toroidal_angles,
+    poloidal_angles,
+    wall_s,
+    radial_build_dict,
+    # Set num_ribs and num_rib_pts such that four (60/12 - 1) additional
+    # locations are interpolated between each specified location
+    num_ribs=61,
+    num_rib_pts=61,
+)
+# Export in-vessel component files
+stellarator.export_invessel_build()
+
+# Construct magnets
+stellarator.construct_magnets(
+    coils_file, width, thickness, toroidal_extent, sample_mod=6
+)
+# Export magnet files
+stellarator.export_magnets()
+
+# Define source mesh parameters
+mesh_size = (11, 81, 61)
+toroidal_extent = 90.0
+# Construct source
+stellarator.construct_source_mesh(mesh_size, toroidal_extent)
+# Export source file
+stellarator.export_source_mesh(filename="source_mesh", export_dir=export_dir)
+
+# Build Cubit model of Parastell Components
+stellarator.build_cubit_model(skip_imprint=False, legacy_faceting=True)
+
+# Export DAGMC neutronics H5M file
+stellarator.export_dagmc(filename="dagmc", export_dir=export_dir)

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(

parastell/magnet_coils.py

@@ -8,7 +8,7 @@
 
 from . import log
 from . import cubit_io as cubit_io
-from .utils import read_yaml_config, filter_kwargs, m2cm
+from .utils import read_yaml_config, filter_kwargs, reorder_loop, m2cm
 
 export_allowed_kwargs = ["step_filename", "export_mesh", "mesh_filename"]
 
@@ -191,18 +191,10 @@ def _filter_coils(self):
             for coil in self.magnet_coils
             if coil.in_toroidal_extent(lower_bound, upper_bound)
         ]
-
-        # Compute toroidal angles of filament centers of mass
-        com_list = np.array([coil.center_of_mass for coil in filtered_coils])
-        com_toroidal_angles = np.arctan2(com_list[:, 1], com_list[:, 0])
-        # Ensure angles are positive
-        com_toroidal_angles = (com_toroidal_angles + 2 * np.pi) % (2 * np.pi)
+        self.magnet_coils = filtered_coils
 
         # Sort coils by center-of-mass toroidal angle and overwrite stored list
-        self.magnet_coils = [
-            coil
-            for _, coil in sorted(zip(com_toroidal_angles, filtered_coils))
-        ]
+        self.magnet_coils = self.sort_coils_toroidally()
 
     def _cut_magnets(self):
         """Cuts the magnets at the planes defining the toriodal extent.
@@ -299,6 +291,17 @@ def export_mesh(self, mesh_filename="magnet_mesh", export_dir=""):
             filename=mesh_filename, export_dir=export_dir
         )
 
+    def sort_coils_toroidally(self):
+        """Reorders list of coils by toroidal angle on range [-pi, pi].
+
+        Arguments:
+            magnet_coils (list of object): list of MagnetCoil class objects.
+
+        Returns:
+            (list of object): sorted list of MagnetCoil class objects.
+        """
+        return sorted(self.magnet_coils, key=lambda x: x.com_toroidal_angle())
+
 
 class MagnetCoil(object):
     """An object representing a single modular stellarator magnet coil.
@@ -340,36 +343,6 @@ def coords(self, data):
             tangents / np.linalg.norm(tangents, axis=1)[:, np.newaxis]
         )
 
-    def in_toroidal_extent(self, lower_bound, upper_bound):
-        """Determines if the coil lies within a given toroidal angular extent,
-        based on filament coordinates.
-
-        Arguments:
-            lower_bound (float): lower bound of toroidal extent [rad].
-            upper_bound (float): upper bound of toroidal extent [rad].
-
-        Returns:
-            in_toroidal_extent (bool): flag to indicate whether coil lies
-                within toroidal bounds.
-        """
-        # Compute toroidal angle of each point in filament
-        toroidal_angles = np.arctan2(self._coords[:, 1], self._coords[:, 0])
-        # Ensure angles are positive
-        toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi)
-        # Compute bounds of toroidal extent of filament
-        min_tor_ang = np.min(toroidal_angles)
-        max_tor_ang = np.max(toroidal_angles)
-
-        # Determine if filament toroidal extent overlaps with that of model
-        if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or (
-            max_tor_ang >= lower_bound or max_tor_ang <= upper_bound
-        ):
-            in_toroidal_extent = True
-        else:
-            in_toroidal_extent = False
-
-        return in_toroidal_extent
-
     def create_magnet(self):
         """Creates a single magnet coil CAD solid in CadQuery.
 
@@ -436,6 +409,83 @@ def create_magnet(self):
         shell = cq.Shell.makeShell(face_list)
         self.solid = cq.Solid.makeSolid(shell)
 
+    def in_toroidal_extent(self, lower_bound, upper_bound):
+        """Determines if the coil lies within a given toroidal angular extent,
+        based on filament coordinates.
+
+        Arguments:
+            lower_bound (float): lower bound of toroidal extent [rad].
+            upper_bound (float): upper bound of toroidal extent [rad].
+
+        Returns:
+            in_toroidal_extent (bool): flag to indicate whether coil lies
+                within toroidal bounds.
+        """
+        # Compute toroidal angle of each point in filament
+        toroidal_angles = np.arctan2(self._coords[:, 1], self._coords[:, 0])
+        # Ensure angles are positive
+        toroidal_angles = (toroidal_angles + 2 * np.pi) % (2 * np.pi)
+        # Compute bounds of toroidal extent of filament
+        min_tor_ang = np.min(toroidal_angles)
+        max_tor_ang = np.max(toroidal_angles)
+
+        # Determine if filament toroidal extent overlaps with that of model
+        if (min_tor_ang >= lower_bound or min_tor_ang <= upper_bound) or (
+            max_tor_ang >= lower_bound or max_tor_ang <= upper_bound
+        ):
+            in_toroidal_extent = True
+        else:
+            in_toroidal_extent = False
+
+        return in_toroidal_extent
+
+    def com_toroidal_angle(self):
+        """Computes the toroidal angle of the coil center of mass, based on
+        filament coordinates.
+
+        Returns:
+            (float): toroidal angle of coil center of mass [rad].
+        """
+        return np.arctan2(self.center_of_mass[1], self.center_of_mass[0])
+
+    def get_ob_mp_index(self):
+        """Finds the index of the outboard midplane coordinate on a coil
+        filament.
+
+        Returns:
+            outboard_index (int): index of the outboard midplane point.
+        """
+        # Compute radial distance of coordinates from z-axis
+        radii = np.linalg.norm(self.coords[:, :2], axis=1)
+        # Determine whether adjacent points cross the midplane (if so, they will
+        # have opposite signs)
+        shifted_coords = np.append(self.coords[1:], [self.coords[1]], axis=0)
+        midplane_flags = -np.sign(self.coords[:, 2] * shifted_coords[:, 2])
+        # Find index of outboard midplane point
+        outboard_index = np.argmax(midplane_flags * radii)
+
+        return outboard_index
+
+    def reorder_coords(self, index):
+        """Reorders coil filament coordinate loop about a given index.
+
+        Arguments:
+            index (int): index about which to reorder coordinate loop.
+        """
+        self.coords = reorder_loop(self.coords, index)
+
+    def orient_coords(self, positive=True):
+        """Orients coil filament coordinate loop such that they initially
+        progress positively or negatively.
+
+        Arguments:
+            positive (bool): progress coordinates in positive direciton
+                (defaults to True). If negative, coordinates will progress in
+                negative direction.
+        """
+        if positive == (self.coords[0, 2] > self.coords[1, 2]):
+            self.coords = np.flip(self.coords, axis=0)
+
 
 def parse_args():
     """Parser for running as a script"""