diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml new file mode 100644 index 0000000..744fe69 --- /dev/null +++ b/.github/workflows/tests.yml @@ -0,0 +1,31 @@ +name: Tests + +on: + push: + branches: + - main + pull_request: + +jobs: + pytest: + runs-on: ubuntu-latest + + strategy: + matrix: + python-version: ["3.10", "3.11", "3.12", "3.13", "3.14"] + + steps: + - name: Checkout + uses: actions/checkout@v4 + + - name: Set up Python + uses: actions/setup-python@v5 + with: + python-version: ${{ matrix.python-version }} + cache: pip + + - name: Install package + run: pip install -e ".[dev]" + + - name: Run tests + run: pytest diff --git a/machine_design/__init__.py b/machine_design/__init__.py index f768f76..e69de29 100644 --- a/machine_design/__init__.py +++ b/machine_design/__init__.py @@ -1,12 +0,0 @@ -from .design import Design -from .generators import ( - FourStupid, - HacklGenerator_3BrokenLines, - HacklGenerator_OneLambda, - HacklGenerator_SixLambdas, - HacklGenerator_TwoLambdas, - save_params, -) -from .geometry import analyze_results, plot_barriers, rotate -from .load import load_design -from .optimization import init_points, objective, objective_single, objective_transform diff --git a/machine_design/designs/__init__.py b/machine_design/designs/__init__.py new file mode 100644 index 0000000..3686109 --- /dev/null +++ b/machine_design/designs/__init__.py @@ -0,0 +1,2 @@ +from .design import Design +from .load import load_design diff --git a/machine_design/designs/computation.py b/machine_design/designs/computation.py new file mode 100644 index 0000000..f75c095 --- /dev/null +++ b/machine_design/designs/computation.py @@ -0,0 +1,113 @@ +from abc import ABC, abstractmethod + +from ansys.aedt.core import Maxwell2d + +from .geometry import GeometryBase + + +class ComputationBase(ABC): + def __init__(self, geometry: GeometryBase) -> None: + self.geometry = geometry + self.setup_name = "Setup1" + self.set_oper_params() + self.set_solution_expressions() + self.set_output_vars() + self.set_post_params() + + @abstractmethod + def set_oper_params(self): ... + + @abstractmethod + def set_solution_expressions(self): ... + + @abstractmethod + def set_output_vars(self): ... + + @abstractmethod + def set_post_params(self): ... + + @abstractmethod + def assign_stator_coils(self, m2d: Maxwell2d) -> None: ... + + @abstractmethod + def inductance_computation(self, m2d: Maxwell2d) -> None: ... + + @abstractmethod + def set_variables(self, m2d: Maxwell2d, *args): ... + + @abstractmethod + def extract_results(self, solutions): ... + + def push_variables(self, m2d: Maxwell2d) -> None: + for k, v in self.oper_params.items(): + m2d[k] = v + + def create_setup(self, m2d: Maxwell2d) -> None: + self.assign_stator_coils(m2d) + self.inductance_computation(m2d) + + # model depth + m2d.model_depth = "StackLength" + # symmetry + m2d.change_symmetry_multiplier("SymmetryFactor") + # Calculation setup + setup = m2d.create_setup(name=self.setup_name) + setup.props["StopTime"] = "Nper/f" + setup.props["TimeStep"] = "1/(f*(PointPer-1))" + setup.props["SaveFieldsType"] = "None" + setup.props["OutputPerObjectCoreLoss"] = False + setup.props["OutputPerObjectSolidLoss"] = True + setup.props["OutputError"] = True + setup.update() + m2d.validate_simple() + + for k, v in self.output_vars.items(): + m2d.create_output_variable(k, v) + + for k, v in self.post_params.items(): + expressions = list(k) if isinstance(k, tuple) else [k] # if multiple report, use list(k). Else, use k + m2d.post.create_report( + expressions=expressions, + setup_sweep_name="", + domain="Sweep", + variations=None, + primary_sweep_variable="Time", + secondary_sweep_variable=None, + report_category=None, + plot_type="Rectangular Plot", + context=None, + subdesign_id=None, + polyline_points=1001, + plot_name=v, + ) + + # TODO: change the other arguments to kwargs + def compute(self, m2d: Maxwell2d, rotor_id, *args, NUM_CORES: int = 1): + assert m2d.mesh is not None + assert m2d.post is not None + + m2d.mesh.assign_length_mesh( + assignment=rotor_id, + inside_selection=True, + maximum_length=3, + maximum_elements=None, + name="rotor", + ) + # core loss rotor + m2d.set_core_losses("Rotor", core_loss_on_field=False) + + self.set_variables(m2d, *args) + + # Analyze + m2d.analyze_setup(self.setup_name, use_auto_settings=False, cores=NUM_CORES) + + solutions = m2d.post.get_solution_data(expressions=self.solution_expressions, primary_sweep_variable="Time") + try: + result = self.extract_results(solutions) + except AttributeError: + result = None + + # Delete solution data to prevent the saving size to explode + m2d.odesign.DeleteFullVariation("All", False) + + return result diff --git a/machine_design/designs/design.py b/machine_design/designs/design.py new file mode 100644 index 0000000..3ec1bff --- /dev/null +++ b/machine_design/designs/design.py @@ -0,0 +1,61 @@ +from ansys.aedt.core import Desktop, Maxwell2d + +from .computation import ComputationBase +from .geometry import GeometryBase + + +class Design: + def __init__(self, m2d: Maxwell2d, geometry: GeometryBase, computation: ComputationBase) -> None: + assert computation.geometry is geometry + self.m2d = m2d + self.geometry = geometry + self.computation = computation + + @classmethod + def create(cls, project_name: str, design_name: str, file_name: str, geometry: GeometryBase, computation: ComputationBase, **kwargs) -> "Design": + m2d = Maxwell2d(project=project_name, design=design_name, solution_type="TransientXY", **kwargs) + obj = cls(m2d, geometry, computation) + obj.add_stator() + obj.save_project(file_name) + return obj + + @classmethod + def load(cls, file_name: str, geometry: GeometryBase, computation: ComputationBase, **kwargs) -> "Design": + desktop = Desktop(**kwargs) + desktop.load_project(file_name) + + m2d = Maxwell2d() + m2d.set_active_design("Design01") + + return cls(m2d, geometry, computation) + + def add_stator(self) -> None: + self.computation.push_variables(self.m2d) + self.geometry.build_stator(self.m2d) + self.computation.create_setup(self.m2d) + + def add_rotor(self) -> None: + self.geometry.build_rotor(self.m2d) + + def add_rotor_barrier(self, barrier_points, segment_type=None) -> None: + self.geometry.add_rotor_barrier(self.m2d, barrier_points, segment_type) + + def delete_rotor(self) -> None: + self.geometry.delete_rotor(self.m2d) + + def compute(self, *args, NUM_CORES: int = 1): + return self.computation.compute(self.m2d, self.geometry.rotor_id, *args, NUM_CORES=NUM_CORES) + + def save_design(self, file_name: str, **kwargs) -> None: + show = kwargs.pop("show", False) + view = kwargs.pop("view", "xy") + self.m2d.plot(show=show, output_file=file_name, view=view) + + def save_project(self, file_name: str | None = None) -> None: + if file_name is None: + self.m2d.save_project() + else: + self.m2d.save_project(file_name) + + def close_project(self) -> None: + self.m2d.close_desktop() diff --git a/machine_design/designs/geometry.py b/machine_design/designs/geometry.py new file mode 100644 index 0000000..58d9dc2 --- /dev/null +++ b/machine_design/designs/geometry.py @@ -0,0 +1,192 @@ +from abc import ABC, abstractmethod + +import numpy as np +from ansys.aedt.core import Maxwell2d +from ansys.aedt.core.modeler.modeler_2d import Modeler2D + + +class GeometryBase(ABC): + def __init__(self) -> None: + self.set_iron() + self.set_geom_params() + self.set_slot_params() + self.set_winds_params() + self.set_mod_params() + self.set_rot_points() + self.set_derived_params() + self.set_udp_par_list_stator() + + @abstractmethod + def set_iron(self): ... + + @abstractmethod + def set_geom_params(self): ... + + @abstractmethod + def set_slot_params(self): ... + + @abstractmethod + def set_winds_params(self): ... + + @abstractmethod + def set_mod_params(self): ... + + @abstractmethod + def set_rot_points(self): ... + + @abstractmethod + def set_derived_params(self): ... + + @abstractmethod + def set_udp_par_list_stator(self): ... + + @abstractmethod + def build_stator(self, m2d: Maxwell2d) -> None: ... + + def mm_to_str(self, var, field) -> float: + val = getattr(self, var)[field] + if not val.endswith("mm"): + raise Exception("val must end with mm") + return float(val[:-2]) + + def _set_appearance(self, obj, color, transparency: float) -> None: + obj.color = color + obj.transparency = transparency + + def set_appearances(self, m2d: Maxwell2d) -> None: + for item in [self.shaft_id, self.region_id, self.band_id]: + self._set_appearance(item, (0, 255, 255), 0.95) + self._set_appearance(self.stator_id, (192, 192, 192), 0.0) + for name in self.id_coils: + self._set_appearance(m2d.modeler[name], (255, 128, 0), 0.0) + self._set_appearance(self.rotor_id, (192, 192, 192), 0.0) + + def build_rotor(self, m2d: Maxwell2d) -> None: + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + rotor_id = modeler.create_polyline(points=self.rot_points, segment_type=["Arc", "Line", "Arc"], cover_surface=True, name="Rotor") + self.rotor_id = rotor_id + rotor_id.material_name = self.Fe + + def add_rotor_barrier(self, m2d: Maxwell2d, barrier_points, segment_type=None) -> None: + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + # Round it for Ansys + barrier_points = np.round(barrier_points, 6) + + # Potentially add the z axis + if barrier_points.shape[1] == 2: + barrier_points = np.hstack((barrier_points, np.zeros((len(barrier_points), 1)))) + + # Convert them into a string format and interpolate + points_str = [[str(y) for y in x] for x in barrier_points] + barrier_id = modeler.create_polyline(points=points_str, segment_type=segment_type, cover_surface=True, name="Barrier") + + # Remove the barrier + self.rotor_id.subtract(barrier_id) + modeler.delete(barrier_id) + + def delete_rotor(self, m2d: Maxwell2d) -> None: + assert isinstance(m2d.modeler, Modeler2D) + m2d.modeler.delete(self.rotor_id) + + def _push_stator_variables(self, m2d: Maxwell2d) -> None: + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + modeler.model_units = "mm" + for k, v in self.geom_params.items(): + m2d[k] = v + for k, v in self.wind_params.items(): + m2d[k] = v + for k, v in self.slot_params.items(): + m2d[k] = v + for k, v in self.mod_params.items(): + m2d[k] = v + + def _build_vacuum_regions(self, m2d: Maxwell2d): + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + # Outer region + region_id = modeler.create_circle( + origin=[0, 0, 0], + radius="DiaStatorYoke/2", + # num_sides="SegAngle", + is_covered=True, + name="Region", + ) + # Band + band_id = modeler.create_circle( + origin=[0, 0, 0], + radius="(DiaStatorGap - (1.0 * Airgap))/2", + # num_sides="mapping_angle", + is_covered=True, + name="Band", + ) + # Shaft + shaft_id = modeler.create_circle( + origin=[0, 0, 0], + radius="DiaShaft/2", + is_covered=True, + name="Shaft", + ) + + # Fit all view + modeler.fit_all() + + return shaft_id, region_id, band_id + + def _split_for_symmetry(self, m2d: Maxwell2d, object_list) -> None: + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + modeler.create_coordinate_system( + origin=[0, 0, 0], + reference_cs="Global", + name="Section", + mode="axis", + x_pointing=["cos(360deg/SymmetryFactor)", "sin(360deg/SymmetryFactor)", 0], + y_pointing=["-sin(360deg/SymmetryFactor)", "cos(360deg/SymmetryFactor)", 0], + ) + modeler.set_working_coordinate_system("Section") + modeler.split(assignment=object_list, plane="ZX", sides="NegativeOnly") + modeler.set_working_coordinate_system("Global") + modeler.split(assignment=object_list, plane="ZX", sides="PositiveOnly") + + def _assign_boundary_conditions(self, m2d: Maxwell2d) -> None: + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + # Symmetrical boundary conditions + pos_1 = "DiaStatorGap/4" + id_bc_1 = modeler.get_edgeid_from_position(position=[pos_1, 0, 0], assignment="Region") + id_bc_2 = modeler.get_edgeid_from_position( + position=[ + pos_1 + "*cos((360deg/SymmetryFactor))", + pos_1 + "*sin((360deg/SymmetryFactor))", + 0, + ], + assignment="Region", + ) + m2d.assign_master_slave( + independent=id_bc_1, + dependent=id_bc_2, + reverse_master=True, + reverse_slave=False, + same_as_master=False, + boundary="Symmetry", + ) + # Zero vector potential + pos_2 = "(DiaStatorYoke/2)" + id_bc_az = modeler.get_edgeid_from_position( + position=[ + pos_2 + "*cos((360deg/SymmetryFactor/2))", + pos_2 + "*sin((360deg/SymmetryFactor)/2)", + 0, + ], + assignment="Region", + ) + m2d.assign_vector_potential(assignment=id_bc_az, vector_value=0, boundary="A0") diff --git a/machine_design/load.py b/machine_design/designs/load.py similarity index 77% rename from machine_design/load.py rename to machine_design/designs/load.py index cd12b03..a6d01a0 100644 --- a/machine_design/load.py +++ b/machine_design/designs/load.py @@ -1,6 +1,8 @@ import os +from .computation import ComputationBase from .design import Design +from .geometry import GeometryBase def load_design( @@ -8,6 +10,8 @@ def load_design( project_name, design_name, aedt_version, + geometry: GeometryBase, + computation: ComputationBase, new_desktop=False, non_graphical=True, close_on_exit=True, @@ -19,6 +23,8 @@ def load_design( project_name, design_name, file_name_aedt, + geometry, + computation, version=aedt_version, non_graphical=non_graphical, new_desktop=new_desktop, @@ -28,6 +34,8 @@ def load_design( else: return Design.load( file_name_aedt, + geometry, + computation, version=aedt_version, non_graphical=non_graphical, new_desktop=new_desktop, diff --git a/machine_design/optimization/__init__.py b/machine_design/optimization/__init__.py new file mode 100644 index 0000000..200594f --- /dev/null +++ b/machine_design/optimization/__init__.py @@ -0,0 +1,10 @@ +from .generators import ( + FourStupid, + HacklGenerator_3BrokenLines, + HacklGenerator_OneLambda, + HacklGenerator_SixLambdas, + HacklGenerator_TwoLambdas, + save_params, +) +from .geometry import analyze_results, plot_barriers, rotate +from .optimization import init_points, objective, objective_single, objective_transform diff --git a/machine_design/generators.py b/machine_design/optimization/generators.py similarity index 99% rename from machine_design/generators.py rename to machine_design/optimization/generators.py index 0798ccb..52612c0 100644 --- a/machine_design/generators.py +++ b/machine_design/optimization/generators.py @@ -6,7 +6,7 @@ from scipy.interpolate import CubicSpline from shapely.geometry import LineString -from .design import Design +from ..designs.design import Design from .geometry import rotate @@ -36,8 +36,8 @@ def save_params(params, file_name: str) -> None: class BarrierGenerator(ABC): # TODO: rename offset def __init__(self, design: Design, r_stator_end: float, offset: None | float = None, n_curve: int = 500, n_flat: int = 20) -> None: - self.r_max = design.rotor_r_max - self.r_min = design.rotor_r_min + self.r_max = design.geometry.rotor_r_max + self.r_min = design.geometry.rotor_r_min self.R = self.r_max - r_stator_end self.offset = offset self.n_curve = n_curve diff --git a/machine_design/geometry.py b/machine_design/optimization/geometry.py similarity index 93% rename from machine_design/geometry.py rename to machine_design/optimization/geometry.py index ecac422..e628a16 100644 --- a/machine_design/geometry.py +++ b/machine_design/optimization/geometry.py @@ -17,7 +17,7 @@ def rotate(x, y, alpha, **kwargs): def plot_barriers(barriers, design, title=None, file_name=None): plt.figure() alphas = np.linspace(0 * np.pi, 2 * np.pi / 4, 100) - for r in [design.rotor_r_min, design.rotor_r_max]: + for r in [design.geometry.rotor_r_min, design.geometry.rotor_r_max]: plt.plot(r * np.cos(alphas), r * np.sin(alphas)) for barrier in barriers: plt.plot(barrier[:, 0], barrier[:, 1]) diff --git a/machine_design/optimization.py b/machine_design/optimization/optimization.py similarity index 100% rename from machine_design/optimization.py rename to machine_design/optimization/optimization.py diff --git a/motors/__init__.py b/motors/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/motors/motor1.py b/motors/motor1.py new file mode 100644 index 0000000..ce1c118 --- /dev/null +++ b/motors/motor1.py @@ -0,0 +1,266 @@ +from ansys.aedt.core import Maxwell2d +from ansys.aedt.core.modeler.modeler_2d import Modeler2D + +from machine_design.designs.computation import ComputationBase +from machine_design.designs.geometry import GeometryBase + + +class Geometry(GeometryBase): + def set_iron(self): + self.Fe = "Cogent Power - M350-50A, B-H at 50Hz" + + def set_geom_params(self): + self.geom_params = { + "DiaStatorGap": "79mm", + "DiaStatorYoke": "125mm", + "Airgap": "0.225mm", + "SlotNumber": "36", + "SlotType": "3", + "DiaShaft": "25mm", + "StackLength": "85mm", + } + + def set_slot_params(self): + self.slot_params = { + "Hs0": "0.95mm", + "Hs1": "0.31mm", + "Hs2": "8.24mm", + "Bs0": "2.2mm", + "Bs1": "3.37mm", + "Bs2": "4.8mm", + "Rs": "1.5mm", + "SetAngle": "10deg", + } + + def set_winds_params(self): + self.wind_params = { + "Layers": "1", + "ParallelPaths": "1", + "CoilPitch": "9", # coil pitch in slots + "SlotLiner": "0.3mm", + "SpaceLayers": "0.2mm", + "Nc": "68", # turns per coil + } + + def set_mod_params(self): + self.PolePairs = 2 + self.mod_params = { + "Poles": f"2*{self.PolePairs}", + "ModelLength": "85mm", + "SymmetryFactor": "Poles", + "StatorSkewAngle": "0deg", + } + + def set_rot_points(self): + self.rot_points = [ + ["DiaShaft/2*cos(360deg/SymmetryFactor)", "DiaShaft/2*sin(360deg/SymmetryFactor)", "0mm"], + ["DiaShaft/2*cos(360deg/(2*SymmetryFactor))", "DiaShaft/2*sin(360deg/(2*SymmetryFactor))", "0mm"], + ["DiaShaft/2", "0mm", "0mm"], + ["DiaStatorGap/2-Airgap", "0mm", "0mm"], + ["(DiaStatorGap/2-Airgap)*cos(360deg/(2*SymmetryFactor))", "(DiaStatorGap/2-Airgap)*sin(360deg/(2*SymmetryFactor))", "0mm"], + ["(DiaStatorGap/2-Airgap)*cos(360deg/SymmetryFactor)", "(DiaStatorGap/2-Airgap)*sin(360deg/SymmetryFactor)", "0mm"], + ] + + def set_derived_params(self): + self.rotor_r_min = self.mm_to_str("geom_params", "DiaShaft") / 2 + self.rotor_r_max = self.mm_to_str("geom_params", "DiaStatorGap") / 2 - self.mm_to_str("geom_params", "Airgap") + + def set_udp_par_list_stator(self): + self.udp_par_list_stator = [ + ["DiaGap", "DiaStatorGap"], + ["DiaYoke", "DiaStatorYoke"], + ["Length", "0mm"], + ["Skew", "0deg"], + ["Slots", "SlotNumber"], + ["SlotType", "SlotType"], + ["Hs0", "Hs0"], + ["Hs01", "0mm"], + ["Hs1", "Hs1"], + ["Hs2", "Hs2"], + ["Bs0", "Bs0"], + ["Bs1", "Bs1"], + ["Bs2", "Bs2"], + ["Rs", "Rs"], + ["FilletType", "0"], + ["HalfSlot", "0"], + ["SegAngle", "0deg"], + ["LenRegion", "0mm"], + ["InfoCore", "0"], + ] + + def build_stator(self, m2d: Maxwell2d) -> None: + self._push_stator_variables(m2d) + shaft_id, region_id, band_id = self._build_vacuum_regions(m2d) + self._assign_rotor_motion(m2d) + stator_id = self._build_stator_core(m2d) + id_coils = self._build_stator_coils(m2d) + self._split_for_symmetry(m2d, [stator_id, shaft_id, region_id, band_id]) + self._assign_boundary_conditions(m2d) + self._assign_stator_mesh(m2d, stator_id, id_coils) + + # core loss + m2d.set_core_losses("Stator", core_loss_on_field=False) + + self.stator_id = stator_id + self.shaft_id = shaft_id + self.region_id = region_id + self.band_id = band_id + self.id_coils = id_coils + + def _assign_rotor_motion(self, m2d: Maxwell2d) -> None: + m2d.assign_rotate_motion( + assignment="Band", + coordinate_system="Global", + axis="Z", + positive_movement=True, + start_position="InitPos", + angular_velocity="RotSpeed", + has_rotation_limits=False, + ) + + def _build_stator_core(self, m2d: Maxwell2d): + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + stator_id = modeler.create_udp( + dll="RMxprt/SlotCore.dll", + parameters=self.udp_par_list_stator, + library="syslib", + name="Stator", + # SolveInside="True", + ) + stator_id.material_name = self.Fe + + return stator_id + + def _build_stator_coils(self, m2d: Maxwell2d): + modeler = m2d.modeler + assert isinstance(modeler, Modeler2D) + + coil_id = modeler.create_rectangle( + origin=["DiaStatorGap/2+Hs0+Hs1+SlotLiner", "-(Bs1/2-SlotLiner)", 0], + sizes=["(Hs2+Rs/2-2*SlotLiner)/Layers-(Layers-1)*SpaceLayers/2", "Bs1-2*SlotLiner", 0], + name="Coil", + material="copper", + ) + modeler.rotate(assignment=coil_id, axis="Z", angle="360deg/SlotNumber/2") + coil_id.duplicate_around_axis(axis="Z", angle="360deg/SlotNumber", clones="SlotNumber/Poles", create_new_objects=True) + id_coils = modeler.get_objects_w_string(string_name="Coil", case_sensitive=True) + + return id_coils + + def _assign_stator_mesh(self, m2d: Maxwell2d, stator_id, id_coils) -> None: + m2d.mesh.assign_length_mesh( + assignment=id_coils, + inside_selection=True, + maximum_length=3, + maximum_elements=None, + name="coils", + ) + m2d.mesh.assign_length_mesh( + assignment=stator_id, + inside_selection=True, + maximum_length=3, + maximum_elements=None, + name="stator", + ) + + +class Computation(ComputationBase): + def set_oper_params(self): + f = 50 # [Hz] + RotSpeed = 60 * f / self.geometry.PolePairs # [rpm] + self.oper_params = { + "Im": "1.5*sqrt(2)A", + "epsI": "pi/4", # current angle + "InitPos": "-30deg", + "f": f"{f}Hz", + "RotSpeed": f"{RotSpeed}rpm", + "Nper": "1/6", # number of included periods + "PointPer": "101", # number of time points per period + } + + def set_solution_expressions(self): + self.solution_expressions = "Moving1.Torque" + + def set_output_vars(self): + self.output_vars = { + "pos": "(Moving1.Position -InitPos) * Poles/2", + "cos0": "cos(pos)", + "cos1": "cos(pos-2*PI/3)", + "cos2": "cos(pos-4*PI/3)", + "sin0": "sin(pos)", + "sin1": "sin(pos-2*PI/3)", + "sin2": "sin(pos-4*PI/3)", + "Lad": "L(PhaseA,PhaseA)*cos0 + L(PhaseA,PhaseB)*cos1 + L(PhaseA,PhaseC)*cos2", + "Laq": "L(PhaseA,PhaseA)*sin0 + L(PhaseA,PhaseB)*sin1 + L(PhaseA,PhaseC)*sin2", + "Lbd": "L(PhaseB,PhaseA)*cos0 + L(PhaseB,PhaseB)*cos1 + L(PhaseB,PhaseC)*cos2", + "Lbq": "L(PhaseB,PhaseA)*sin0 + L(PhaseB,PhaseB)*sin1 + L(PhaseB,PhaseC)*sin2", + "Lcd": "L(PhaseC,PhaseA)*cos0 + L(PhaseC,PhaseB)*cos1 + L(PhaseC,PhaseC)*cos2", + "Lcq": "L(PhaseC,PhaseA)*sin0 + L(PhaseC,PhaseB)*sin1 + L(PhaseC,PhaseC)*sin2", + "L_d": "(Lad*cos0 + Lbd*cos1 + Lcd*cos2) * 2/3", + "L_q": "(Laq*sin0 + Lbq*sin1 + Lcq*sin2) * 2/3", + "Flux_d": "(FluxLinkage(PhaseA)*cos0+FluxLinkage(PhaseB)*cos1+FluxLinkage(PhaseC)*cos2)*2/3", + "Flux_q": "-(FluxLinkage(PhaseA)*sin0+FluxLinkage(PhaseB)*sin1+FluxLinkage(PhaseC)*sin2)*2/3", + "Ui_d": "(InducedVoltage(PhaseA)*cos0+InducedVoltage(PhaseB)*cos1+InducedVoltage(PhaseC)*cos2)*2/3", + "Ui_q": "-(InducedVoltage(PhaseA)*sin0+InducedVoltage(PhaseB)*sin1+InducedVoltage(PhaseC)*sin2)*2/3", + "I_d": "(InputCurrent(PhaseA)*cos0 + InputCurrent(PhaseB)*cos1 + InputCurrent(PhaseC)*cos2)*2/3", + "I_q": "-(InputCurrent(PhaseA)*sin0 + InputCurrent(PhaseB)*sin1 + InputCurrent(PhaseC)*sin2)*2/3", + "Irms": "sqrt(I_d^2+I_q^2)/sqrt(2)", + } + + def set_post_params(self): + self.post_params = { # reports + ("InducedVoltage(PhaseA)", "InducedVoltage(PhaseB)", "InducedVoltage(PhaseC)"): "InducedVoltage", + ("Moving1.Torque"): "Torque", + ("InputCurrent(PhaseA)", "InputCurrent(PhaseB)", "InputCurrent(PhaseC)"): "Current", + ( + "FluxLinkage(PhaseA)", + "FluxLinkage(PhaseB)", + "FluxLinkage(PhaseC)", + ): "FluxLinkage", + ("I_d", "I_q"): "Current_dq", + ("Flux_d", "Flux_q"): "FluxLinkage_dq", + ("Ui_d", "Ui_q"): "InducedVoltage_dq", + ("L_d", "L_q"): "Inductance_dq", + } + + def assign_stator_coils(self, m2d: Maxwell2d) -> None: + # Excitations + I_A = "Im * cos(2*pi*f*time+epsI)" + I_B = "Im * cos(2*pi*f*time-120deg+epsI)" + I_C = "Im * cos(2*pi*f*time-240deg+epsI)" + + phases_polarity = ["Positive", "Negative", "Positive"] + phases_name = ["A", "C", "B"] + phases_current = [I_A, I_C, I_B] + i_coil = 0 + for phase_polarity, phase_name, phase_current in zip(phases_polarity, phases_name, phases_current): + names = [] + for _ in range(3): + m2d.assign_coil( + assignment=[self.geometry.id_coils[i_coil]], + conductors_number="Nc", + polarity=phase_polarity, + name=f"CS{i_coil + 1}", + ) + names.append(f"CS{i_coil + 1}") + i_coil += 1 + m2d.assign_winding( + assignment=None, + winding_type="Current", + is_solid=False, + current=phase_current, + parallel_branches="ParallelPaths", + name=f"Phase{phase_name}", + ) + m2d.add_winding_coils(assignment=f"Phase{phase_name}", coils=names) + + def inductance_computation(self, m2d: Maxwell2d) -> None: + m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=False) + + def set_variables(self, m2d: Maxwell2d, *args): + pass + + def extract_results(self, solutions): + return solutions.data_magnitude() diff --git a/machine_design/design2.py b/motors/motor2.py similarity index 89% rename from machine_design/design2.py rename to motors/motor2.py index 92fecf4..40c5bd6 100644 --- a/machine_design/design2.py +++ b/motors/motor2.py @@ -1,9 +1,11 @@ -import numpy as np +from ansys.aedt.core import Maxwell2d -from .design import Design +from machine_design.designs.computation import ComputationBase +from .motor1 import Geometry -class Design2(Design): + +class Geometry2(Geometry): def set_geom_params(self): super().set_geom_params() self.geom_params["SlotNumber"] = "40" @@ -18,9 +20,14 @@ def set_winds_params(self): super().set_winds_params() self.wind_params["Nc"] = "113" + def set_derived_params(self): + pass + + +class Computation2(ComputationBase): def set_oper_params(self): f = 50 # [Hz] - RotSpeed = 60 * f / self.PolePairs # [rpm] + RotSpeed = 60 * f / self.geometry.PolePairs # [rpm] self.oper_params = { "Id1": "0.0A", "Iq1": "0.0A", @@ -37,9 +44,6 @@ def set_oper_params(self): "PointPer": "101", # number of time points per period } - def set_derived_params(self): - pass - def set_solution_expressions(self): self.solution_expressions = [ "V_d1", @@ -170,10 +174,28 @@ def set_output_vars(self): def set_post_params(self): self.post_params = { # reports - ("InducedVoltage(PhaseA)", "InducedVoltage(PhaseB)", "InducedVoltage(PhaseC)", "InducedVoltage(PhaseD)", "InducedVoltage(PhaseE)"): "InducedVoltage", + ( + "InducedVoltage(PhaseA)", + "InducedVoltage(PhaseB)", + "InducedVoltage(PhaseC)", + "InducedVoltage(PhaseD)", + "InducedVoltage(PhaseE)", + ): "InducedVoltage", ("Moving1.Torque"): "Torque", - ("InputCurrent(PhaseA)", "InputCurrent(PhaseB)", "InputCurrent(PhaseC)", "InputCurrent(PhaseD)", "InputCurrent(PhaseE)"): "Current", - ("FluxLinkage(PhaseA)", "FluxLinkage(PhaseB)", "FluxLinkage(PhaseC)", "FluxLinkage(PhaseD)", "FluxLinkage(PhaseE)"): "FluxLinkage", + ( + "InputCurrent(PhaseA)", + "InputCurrent(PhaseB)", + "InputCurrent(PhaseC)", + "InputCurrent(PhaseD)", + "InputCurrent(PhaseE)", + ): "Current", + ( + "FluxLinkage(PhaseA)", + "FluxLinkage(PhaseB)", + "FluxLinkage(PhaseC)", + "FluxLinkage(PhaseD)", + "FluxLinkage(PhaseE)", + ): "FluxLinkage", ("I_d1", "I_q1", "I_d3", "I_q3"): "Current_dq", ("Flux_d1", "Flux_q1", "Flux_d3", "Flux_q3"): "FluxLinkage_dq", ("Flux_e_d1", "Flux_e_q1", "Flux_e_d3", "Flux_e_q3"): "FluxLinkage excitation_dq", @@ -183,9 +205,7 @@ def set_post_params(self): ("Ld1q1", "Ld1d3", "Ld1q3", "Lq1d3", "Lq1q3", "Ld3q3"): "Inductance_dq cross-coupling", } - def assign_stator_coils(self): - m2d = self.m2d - + def assign_stator_coils(self, m2d: Maxwell2d) -> None: # Excitations I_A = "Im1*cos(2*pi*f*time+epsI1-pi) + Im3*cos(3*(2*pi*f*time)+epsI3-pi)" I_B = "Im1*cos(2*pi*f*time-72deg+epsI1-pi) + Im3*cos(3*(2*pi*f*time-72deg)+epsI3-pi)" @@ -302,24 +322,25 @@ def assign_stator_coils(self): m2d.add_winding_coils(assignment="PhaseD", coils=["CS2", "CS3"]) m2d.add_winding_coils(assignment="PhaseE", coils=["CS6", "CS7"]) - def inductance_computation(self): - self.m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=True) + def inductance_computation(self, m2d: Maxwell2d) -> None: + m2d.change_inductance_computation(compute_transient_inductance=True, incremental_matrix=True) - def set_variables(self, Id1, Iq1, Id3, Iq3): - self.m2d.variable_manager["Id1"] = f"{Id1}A" - self.m2d.variable_manager["Iq1"] = f"{Iq1}A" - self.m2d.variable_manager["Id3"] = f"{Id3}A" - self.m2d.variable_manager["Iq3"] = f"{Iq3}A" + def set_variables(self, m2d: Maxwell2d, Id1, Iq1, Id3, Iq3): + m2d.variable_manager["Id1"] = f"{Id1}A" + m2d.variable_manager["Iq1"] = f"{Iq1}A" + m2d.variable_manager["Id3"] = f"{Id3}A" + m2d.variable_manager["Iq3"] = f"{Iq3}A" def extract_results(self, solutions): - # TODO: this works only because torque is the last one in the array - out = np.zeros(len(self.solution_expressions)) - for i, expr in enumerate(self.solution_expressions): - data = solutions.data_real(expr) - val = float(np.mean(data[:-1])) + out = {} + for expr in self.solution_expressions: + val = solutions.data_real(expr) - if expr.startswith("L_"): + if expr.startswith("Ld") or expr.startswith("Lq"): val /= 1e9 + elif expr.startswith("I_"): + val /= 1e3 - out[i] = val - return data + # TODO: possibly assign val[:-1]. check whether values are identical + out[expr] = val + return out diff --git a/notebooks/run.py b/notebooks/run.py index 1686e05..34728fe 100644 --- a/notebooks/run.py +++ b/notebooks/run.py @@ -5,15 +5,16 @@ import numpy as np import pandas as pd -from machine_design import ( +from machine_design.designs import load_design +from machine_design.optimization import ( HacklGenerator_3BrokenLines, HacklGenerator_OneLambda, HacklGenerator_SixLambdas, analyze_results, - load_design, plot_barriers, save_params, ) +from motors.motor1 import Computation, Geometry aedt_version = "2024.1" n_designs = 50 @@ -30,7 +31,9 @@ os.makedirs(path, exist_ok=True) file_name_aedt = f"{path_data}/{project_name}.aedt" -design = load_design(file_name_aedt, project_name, design_name, aedt_version) +geometry = Geometry() +computation = Computation(geometry) +design = load_design(file_name_aedt, project_name, design_name, aedt_version, geometry, computation) generators = [ HacklGenerator_OneLambda(design, r_stator_end, offset=offset), HacklGenerator_SixLambdas(design, r_stator_end, offset=offset), diff --git a/notebooks/run_optimization.py b/notebooks/run_optimization.py index 1e61d1a..b32a971 100644 --- a/notebooks/run_optimization.py +++ b/notebooks/run_optimization.py @@ -13,15 +13,16 @@ from botorch.utils.transforms import normalize, unnormalize from gpytorch.mlls import ExactMarginalLogLikelihood -from machine_design import ( +from machine_design.designs import load_design +from machine_design.optimization import ( HacklGenerator_3BrokenLines, HacklGenerator_OneLambda, HacklGenerator_SixLambdas, init_points, - load_design, objective, objective_transform, ) +from motors.motor1 import Computation, Geometry torch.set_default_dtype(torch.float64) @@ -44,7 +45,9 @@ os.makedirs(path_data, exist_ok=True) file_name_aedt = f"{path_data}/{project_name}.aedt" -design = load_design(file_name_aedt, project_name, design_name, aedt_version) +geometry = Geometry() +computation = Computation(geometry) +design = load_design(file_name_aedt, project_name, design_name, aedt_version, geometry, computation) generators = [ HacklGenerator_OneLambda(design, r_stator_end, offset=offset), HacklGenerator_SixLambdas(design, r_stator_end, offset=offset), diff --git a/pyproject.toml b/pyproject.toml index 9a8a1f5..d6e4e24 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -29,12 +29,30 @@ dependencies=[ 'trame', 'matplotlib', 'botorch', + 'pandas', + 'scipy', + 'shapely', +] + +[project.optional-dependencies] +dev = [ + 'pytest', ] [tool.setuptools] packages = [ "machine_design", + "machine_design.designs", + "machine_design.optimization", + "motors", +] + +[tool.pytest.ini_options] +testpaths = ["tests"] +markers = [ + "ansys: requires a running Ansys Electronics Desktop session and a license seat; excluded by default", ] +addopts = "-m 'not ansys'" [tool.ruff] line-length = 250 diff --git a/machine_design/design.py b/tests/legacy/design.py similarity index 99% rename from machine_design/design.py rename to tests/legacy/design.py index 3324a50..bd0ba5d 100644 --- a/machine_design/design.py +++ b/tests/legacy/design.py @@ -1,3 +1,5 @@ +# Frozen snapshot of machine_design/design.py at commit 086817e. +# Golden reference for tests/test_ansys_regression.py - do not update when design.py changes. import numpy as np from ansys.aedt.core import Desktop, Maxwell2d from ansys.aedt.core.modeler.modeler_2d import Modeler2D diff --git a/tests/test_ansys_regression.py b/tests/test_ansys_regression.py new file mode 100644 index 0000000..0e87a54 --- /dev/null +++ b/tests/test_ansys_regression.py @@ -0,0 +1,112 @@ +"""Golden-master check: tests/legacy/design.py vs. the live machine_design/design.py. + +Generates two feasible rotor designs (different seeds) with HacklGenerator_OneLambda and +runs each through add_rotor -> compute -> delete_rotor on the SAME Design instance, both +with the frozen reference implementation and the current one - mirroring how an +optimization loop reuses one Design across many candidate rotors. Asserts the resulting +torque waveforms match at each step. Run explicitly with: + + pytest -m ansys tests/test_ansys_regression.py + +Requires a running Ansys Electronics Desktop session and a free license seat. +""" + +import importlib.util +from pathlib import Path + +import numpy as np +import pytest + +from machine_design.designs.design import Design as LiveDesign +from machine_design.optimization.generators import HacklGenerator_OneLambda +from motors.motor1 import Computation, Geometry + +pytestmark = pytest.mark.ansys + +AEDT_VERSION = "2024.1" +R_STATOR_END = 0.7 +OFFSET = 0.35 +SEEDS = (42, 43) +NUM_CORES = 1 + + +def _load_legacy_design_class(): + legacy_path = Path(__file__).parent / "legacy" / "design.py" + spec = importlib.util.spec_from_file_location("legacy_design", legacy_path) + module = importlib.util.module_from_spec(spec) + spec.loader.exec_module(module) + return module.Design + + +def _generate_one_lambda_barriers(seed): + # rotor_r_min/rotor_r_max are pure-Python derived params, so a throwaway + # instance (no Ansys) is enough to size the generator. + geometry = Geometry() + dummy = LiveDesign(m2d=None, geometry=geometry, computation=Computation(geometry)) + generator = HacklGenerator_OneLambda(dummy, R_STATOR_END, offset=OFFSET) + + np.random.seed(seed) + while True: + params = generator.random_parameters() + generator.set_parameters(params) + barriers = generator.generate_barriers() + barriers = generator.split_barriers(barriers) + if generator.feasible_barriers(barriers): + return barriers + + +def _add_rotor(design, barriers): + design.add_rotor() + for barrier in barriers: + design.add_rotor_barrier(barrier) + + +def _create_design(design_cls, project_name, file_name, extra_args=()): + return design_cls.create( + project_name, + "Design01", + file_name, + *extra_args, + version=AEDT_VERSION, + non_graphical=True, + new_desktop=False, + close_on_exit=False, + ) + + +def _compute_torque(design, barriers, add_rotor): + add_rotor(design, barriers) + torque = design.compute(NUM_CORES=NUM_CORES) + design.delete_rotor() + return torque + + +def test_legacy_and_live_design_match_across_rotor_rebuilds(tmp_path): + barrier_sets = [_generate_one_lambda_barriers(seed) for seed in SEEDS] + LegacyDesign = _load_legacy_design_class() + + legacy_design, live_design = None, None + try: + legacy_design = _create_design(LegacyDesign, "RegressionTest_legacy", str(tmp_path / "legacy.aedt")) + + live_geometry = Geometry() + live_computation = Computation(live_geometry) + live_design = _create_design( + LiveDesign, + "RegressionTest_live", + str(tmp_path / "live.aedt"), + extra_args=(live_geometry, live_computation), + ) + + for seed, barriers in zip(SEEDS, barrier_sets): + torque_legacy = _compute_torque(legacy_design, barriers, _add_rotor) + torque_live = _compute_torque(live_design, barriers, _add_rotor) + + assert torque_legacy is not None, f"legacy torque is None for seed {seed}" + assert torque_live is not None, f"live torque is None for seed {seed}" + np.testing.assert_allclose(torque_live, torque_legacy, rtol=1e-6, err_msg=f"mismatch for seed {seed}") + finally: + if live_design is not None: + live_design.close_project() + elif legacy_design is not None: + legacy_design.close_project() diff --git a/tests/test_design.py b/tests/test_design.py new file mode 100644 index 0000000..bdeae90 --- /dev/null +++ b/tests/test_design.py @@ -0,0 +1,150 @@ +import re + +import pytest + +from motors.motor1 import Computation, Geometry + + +@pytest.fixture +def geometry(): + return Geometry() + + +@pytest.fixture +def computation(geometry): + return Computation(geometry) + + +def test_setup_name(computation): + assert computation.setup_name == "Setup1" + + +def test_iron(geometry): + assert geometry.Fe == "Cogent Power - M350-50A, B-H at 50Hz" + + +def test_geom_params(geometry): + assert geometry.geom_params == { + "DiaStatorGap": "79mm", + "DiaStatorYoke": "125mm", + "Airgap": "0.225mm", + "SlotNumber": "36", + "SlotType": "3", + "DiaShaft": "25mm", + "StackLength": "85mm", + } + + +def test_slot_params(geometry): + assert geometry.slot_params == { + "Hs0": "0.95mm", + "Hs1": "0.31mm", + "Hs2": "8.24mm", + "Bs0": "2.2mm", + "Bs1": "3.37mm", + "Bs2": "4.8mm", + "Rs": "1.5mm", + "SetAngle": "10deg", + } + + +def test_wind_params(geometry): + assert geometry.wind_params == { + "Layers": "1", + "ParallelPaths": "1", + "CoilPitch": "9", + "SlotLiner": "0.3mm", + "SpaceLayers": "0.2mm", + "Nc": "68", + } + + +def test_mod_params(geometry): + assert geometry.PolePairs == 2 + assert geometry.mod_params == { + "Poles": "2*2", + "ModelLength": "85mm", + "SymmetryFactor": "Poles", + "StatorSkewAngle": "0deg", + } + + +def test_oper_params(computation): + assert computation.oper_params == { + "Im": "1.5*sqrt(2)A", + "epsI": "pi/4", + "InitPos": "-30deg", + "f": "50Hz", + "RotSpeed": "1500.0rpm", + "Nper": "1/6", + "PointPer": "101", + } + + +def test_rot_points(geometry): + assert len(geometry.rot_points) == 6 + assert all(len(point) == 3 for point in geometry.rot_points) + + +def test_derived_params(geometry): + assert geometry.rotor_r_min == pytest.approx(12.5) + assert geometry.rotor_r_max == pytest.approx(39.5 - 0.225) + + +def test_solution_expressions(computation): + assert computation.solution_expressions == "Moving1.Torque" + + +def test_udp_par_list_stator(geometry): + keys = [item[0] for item in geometry.udp_par_list_stator] + assert keys == [ + "DiaGap", + "DiaYoke", + "Length", + "Skew", + "Slots", + "SlotType", + "Hs0", + "Hs01", + "Hs1", + "Hs2", + "Bs0", + "Bs1", + "Bs2", + "Rs", + "FilletType", + "HalfSlot", + "SegAngle", + "LenRegion", + "InfoCore", + ] + + +def test_output_vars_only_reference_earlier_keys(computation): + # Ansys evaluates output variables in insertion order, so a formula referencing + # another output variable must not reference one that is defined later. + defined = set() + all_keys = set(computation.output_vars.keys()) + for key, formula in computation.output_vars.items(): + later_keys = all_keys - defined - {key} + referenced_later = [k for k in later_keys if re.search(rf"\b{re.escape(k)}\b", formula)] + assert not referenced_later, f"'{key}' formula references not-yet-defined {referenced_later}" + defined.add(key) + + +def test_output_vars_known_formula(computation): + assert computation.output_vars["Irms"] == "sqrt(I_d^2+I_q^2)/sqrt(2)" + + +def test_post_params_plot_names_unique(computation): + plot_names = list(computation.post_params.values()) + assert len(plot_names) == len(set(plot_names)) + + +def test_mm_to_str(geometry): + assert geometry.mm_to_str("geom_params", "DiaShaft") == pytest.approx(25.0) + + +def test_mm_to_str_raises_without_mm_suffix(geometry): + with pytest.raises(Exception): + geometry.mm_to_str("geom_params", "SlotNumber") diff --git a/tests/test_design2.py b/tests/test_design2.py new file mode 100644 index 0000000..56d6cff --- /dev/null +++ b/tests/test_design2.py @@ -0,0 +1,79 @@ +import pytest + +from motors.motor1 import Geometry +from motors.motor2 import Computation2, Geometry2 + + +@pytest.fixture +def base_geometry(): + return Geometry() + + +@pytest.fixture +def geometry2(): + return Geometry2() + + +@pytest.fixture +def computation2(geometry2): + return Computation2(geometry2) + + +def test_geom_params_overrides_slot_number_only(base_geometry, geometry2): + expected = dict(base_geometry.geom_params) + expected["SlotNumber"] = "40" + assert geometry2.geom_params == expected + + +def test_slot_params_overrides_bs1_bs2_setangle_only(base_geometry, geometry2): + expected = dict(base_geometry.slot_params) + expected["Bs1"] = "3.0mm" + expected["Bs2"] = "4.3mm" + expected["SetAngle"] = "9deg" + assert geometry2.slot_params == expected + + +def test_wind_params_overrides_nc_only(base_geometry, geometry2): + expected = dict(base_geometry.wind_params) + expected["Nc"] = "113" + assert geometry2.wind_params == expected + + +def test_oper_params_fully_replaced(computation2): + assert computation2.oper_params == { + "Id1": "0.0A", + "Iq1": "0.0A", + "Id3": "0.0A", + "Iq3": "0.0A", + "epsI1": "atan2(Iq1,Id1)", + "epsI3": "atan2(Iq3,Id3)", + "Im1": "sqrt(Id1^2+Iq1^2)", + "Im3": "sqrt(Id3^2+Iq3^2)", + "InitPos": "-45deg", + "f": "50Hz", + "RotSpeed": "1500.0rpm", + "Nper": "1/10", + "PointPer": "101", + } + + +def test_derived_params_are_not_set(geometry2): + # Geometry2.set_derived_params() is a no-op, unlike Geometry's. + assert not hasattr(geometry2, "rotor_r_min") + assert not hasattr(geometry2, "rotor_r_max") + + +def test_solution_expressions_is_a_list_ending_in_torque(computation2): + assert isinstance(computation2.solution_expressions, list) + assert computation2.solution_expressions[-1] == "Moving1.Torque" + assert len(computation2.solution_expressions) == 23 + + +def test_post_params_plot_names_unique(computation2): + plot_names = list(computation2.post_params.values()) + assert len(plot_names) == len(set(plot_names)) + + +def test_output_vars_known_values(computation2): + assert computation2.output_vars["Rstat"] == "19" + assert computation2.output_vars["Lew"] == "0"