diff --git a/machine_design/designs/computation.py b/machine_design/designs/computation.py index 4e57fbe..f756403 100644 --- a/machine_design/designs/computation.py +++ b/machine_design/designs/computation.py @@ -1,5 +1,6 @@ from abc import ABC, abstractmethod +import numpy as np from ansys.aedt.core import Maxwell2d from machine_design.generic.compute_initpos import compute_initpos @@ -43,6 +44,10 @@ def set_variables(self, m2d: Maxwell2d, *args): ... @abstractmethod def extract_results(self, solutions): ... + @staticmethod + def extract_expression(solutions, expr: str) -> np.ndarray: + return np.asarray(solutions.data_real(expr, convert_to_SI=True)) + def set_initpos(self) -> None: Q = int(self.geometry.geom_params["SlotNumber"]) p = self.geometry.PolePairs diff --git a/motors/synrm_3f_36s.py b/motors/synrm_3f_36s.py index 812a051..ff51bbf 100644 --- a/motors/synrm_3f_36s.py +++ b/motors/synrm_3f_36s.py @@ -254,21 +254,20 @@ def set_variables(self, m2d: Maxwell2d, Id, Iq): m2d.variable_manager["Iq"] = f"{Iq}A" def extract_results(self, solutions): - position = np.array(solutions.data_real("Moving1.Position")) - torque = np.array(solutions.data_real("Moving1.Torque")) + # SI units + position = self.extract_expression(solutions, "Moving1.Position") + torque = self.extract_expression(solutions, "Moving1.Torque") - theta_el = np.deg2rad(position - self.InitPos) * self.geometry.PolePairs - flux_phases = np.stack([np.array(solutions.data_real(f"FluxLinkage(Phase{p})")) for p in "ABC"], axis=-1) - vind_phases = np.stack([np.array(solutions.data_real(f"InducedVoltage(Phase{p})")) for p in "ABC"], axis=-1) - current_phases = np.stack([np.array(solutions.data_real(f"InputCurrent(Phase{p})")) for p in "ABC"], axis=-1) + theta_el = (position - np.deg2rad(self.InitPos)) * self.geometry.PolePairs + flux_phases = np.stack([self.extract_expression(solutions, f"FluxLinkage(Phase{p})") for p in "ABC"], axis=-1) + vind_phases = np.stack([self.extract_expression(solutions, f"InducedVoltage(Phase{p})") for p in "ABC"], axis=-1) + current_phases = np.stack([self.extract_expression(solutions, f"InputCurrent(Phase{p})") for p in "ABC"], axis=-1) Flux_d, Flux_q = to_dq(flux_phases, theta_el, harmonic=1) Ui_d, Ui_q = to_dq(vind_phases, theta_el, harmonic=1) - I_d, I_q = (v / 1e3 for v in to_dq(current_phases, theta_el, harmonic=1)) - # mA to A + I_d, I_q = to_dq(current_phases, theta_el, harmonic=1) - L_raw = [np.stack([np.array(solutions.data_real(f"L(Phase{x},Phase{y})")) for y in "ABC"], axis=-1) / 1e9 for x in "ABC"] - # nH to H + L_raw = [np.stack([self.extract_expression(solutions, f"L(Phase{x},Phase{y})") for y in "ABC"], axis=-1) for x in "ABC"] L_d_row = np.zeros((len(position), 3)) L_q_row = np.zeros((len(position), 3)) diff --git a/motors/synrm_3f_60s.py b/motors/synrm_3f_60s.py index 65c0189..ef7b2af 100644 --- a/motors/synrm_3f_60s.py +++ b/motors/synrm_3f_60s.py @@ -1,6 +1,5 @@ """Anchor `synrm_3f_60s`: 3-phase synchronous reluctance machine, 60 stator slots.""" -import numpy as np from ansys.aedt.core import Maxwell2d from machine_design.generic.winding import phase_groups @@ -58,4 +57,4 @@ def assign_stator_coils(self, m2d: Maxwell2d) -> None: m2d.add_winding_coils(assignment=f"Phase{phase_name}", coils=[f"CS_{coil_name}" for coil_name, _ in group]) def extract_results(self, solutions): - return {"Moving1.Torque": np.array(solutions.data_real("Moving1.Torque"))} + return {"Moving1.Torque": self.extract_expression(solutions, "Moving1.Torque")} diff --git a/motors/synrm_5f_40s.py b/motors/synrm_5f_40s.py index d22601a..2242c64 100644 --- a/motors/synrm_5f_40s.py +++ b/motors/synrm_5f_40s.py @@ -223,21 +223,21 @@ def set_variables(self, m2d: Maxwell2d, Id1, Iq1, Id3, Iq3): m2d.variable_manager["Iq3"] = f"{Iq3}A" def extract_results(self, solutions): - position = np.array(solutions.data_real("Moving1.Position")) - torque = np.array(solutions.data_real("Moving1.Torque")) + # SI units + position = self.extract_expression(solutions, "Moving1.Position") + torque = self.extract_expression(solutions, "Moving1.Torque") - theta_el = electrical_angle(position, self.InitPos, self.geometry.PolePairs, self.RotSign, degrees=True) - flux_phases = np.stack([np.array(solutions.data_real(f"FluxLinkage(Phase{p})")) for p in "ABCDE"], axis=-1) - vind_phases = np.stack([np.array(solutions.data_real(f"InducedVoltage(Phase{p})")) for p in "ABCDE"], axis=-1) - current_phases = np.stack([np.array(solutions.data_real(f"InputCurrent(Phase{p})")) for p in "ABCDE"], axis=-1) + theta_el = electrical_angle(position, np.deg2rad(self.InitPos), self.geometry.PolePairs, self.RotSign, degrees=False) + flux_phases = np.stack([self.extract_expression(solutions, f"FluxLinkage(Phase{p})") for p in "ABCDE"], axis=-1) + vind_phases = np.stack([self.extract_expression(solutions, f"InducedVoltage(Phase{p})") for p in "ABCDE"], axis=-1) + current_phases = np.stack([self.extract_expression(solutions, f"InputCurrent(Phase{p})") for p in "ABCDE"], axis=-1) Flux_d1, Flux_q1 = to_dq(flux_phases, theta_el, harmonic=1) Flux_d3, Flux_q3 = to_dq(flux_phases, theta_el, harmonic=3) Vind_d1, Vind_q1 = to_dq(vind_phases, theta_el, harmonic=1) Vind_d3, Vind_q3 = to_dq(vind_phases, theta_el, harmonic=3) - I_d1, I_q1 = (v / 1e3 for v in to_dq(current_phases, theta_el, harmonic=1)) - I_d3, I_q3 = (v / 1e3 for v in to_dq(current_phases, theta_el, harmonic=3)) - # mA to A + I_d1, I_q1 = to_dq(current_phases, theta_el, harmonic=1) + I_d3, I_q3 = to_dq(current_phases, theta_el, harmonic=3) Im1 = np.sqrt(self.Id1**2 + self.Iq1**2) Im3 = np.sqrt(self.Id3**2 + self.Iq3**2) @@ -255,8 +255,7 @@ def extract_results(self, solutions): V_d1, V_q1 = to_dq(V_phases, theta_el, harmonic=1) V_d3, V_q3 = to_dq(V_phases, theta_el, harmonic=3) - L_raw = [np.stack([np.array(solutions.data_real(f"L(Phase{x},Phase{y})")) for y in "ABCDE"], axis=-1) / 1e9 for x in "ABCDE"] - # nH to H + L_raw = [np.stack([self.extract_expression(solutions, f"L(Phase{x},Phase{y})") for y in "ABCDE"], axis=-1) for x in "ABCDE"] L_d1_row = np.zeros((len(time), 5)) L_q1_row = np.zeros((len(time), 5)) diff --git a/motors/synrm_5f_60s.py b/motors/synrm_5f_60s.py index 8e6fe46..9a47438 100644 --- a/motors/synrm_5f_60s.py +++ b/motors/synrm_5f_60s.py @@ -4,7 +4,6 @@ phases. Excitation is dq1 + dq3. """ -import numpy as np from ansys.aedt.core import Maxwell2d from machine_design.generic.winding import phase_groups @@ -64,4 +63,4 @@ def assign_stator_coils(self, m2d: Maxwell2d) -> None: m2d.add_winding_coils(assignment=f"Phase{phase_name}", coils=[f"CS_{coil_name}" for coil_name, _ in group]) def extract_results(self, solutions): - return {"Moving1.Torque": np.array(solutions.data_real("Moving1.Torque"))} + return {"Moving1.Torque": self.extract_expression(solutions, "Moving1.Torque")} diff --git a/tests/test_ansys_regression.py b/tests/test_ansys_regression.py index 5bf8111..6e4c2b4 100644 --- a/tests/test_ansys_regression.py +++ b/tests/test_ansys_regression.py @@ -107,12 +107,7 @@ def test_legacy_and_live_design_match_across_rotor_rebuilds(tmp_path): def _legacy_extract_results(solutions): out = {} for expr in legacy_design.solution_expressions: - val = np.array(solutions.data_real(expr)) - if expr in ("I_d", "I_q"): - val = val / 1e3 # mA to A - elif expr in ("L_d", "L_q"): - val = val / 1e9 # nH to H - out[expr] = val + out[expr] = np.array(solutions.data_real(expr, convert_to_SI=True)) # use SI units to compare return out legacy_design.extract_results = _legacy_extract_results diff --git a/tests/test_synrm_3f_36s.py b/tests/test_synrm_3f_36s.py index 69d9ad5..95ffa38 100644 --- a/tests/test_synrm_3f_36s.py +++ b/tests/test_synrm_3f_36s.py @@ -149,17 +149,18 @@ def test_mm_to_str_raises_without_mm_suffix(geometry): def test_extract_results_converts_units(computation): computation.set_variables(SimpleNamespace(variable_manager={}), Id=1.0, Iq=2.0) - def data_real(expr): + def data_real(expr, convert_to_SI=False): + assert convert_to_SI, f"expected convert_to_SI=True for {expr}" letters = ["ABC".index(part[0]) + 1 for part in expr.split("Phase")[1:]] value = float(np.prod(letters) if letters else 7.0) - return np.full(2, value * 1e9 if expr.startswith("L(") else value) + return np.full(2, value) solutions = SimpleNamespace(data_real=data_real, primary_sweep_values=[0.0, 1.0]) out = computation.extract_results(solutions) - theta_el = np.deg2rad(np.full(2, 7.0) - computation.InitPos) * computation.geometry.PolePairs + theta_el = (np.full(2, 7.0) - np.deg2rad(computation.InitPos)) * computation.geometry.PolePairs current_phases = np.stack([np.full(2, i + 1.0) for i in range(3)], axis=-1) expected_I_d, _ = to_dq(current_phases, theta_el, harmonic=1) - assert out["I_d"] == pytest.approx(expected_I_d / 1e3) + assert out["I_d"] == pytest.approx(expected_I_d) assert out["Moving1.Torque"] == pytest.approx(np.full(2, 7.0)) diff --git a/tests/test_synrm_5f_40s.py b/tests/test_synrm_5f_40s.py index ed7953e..dbfe673 100644 --- a/tests/test_synrm_5f_40s.py +++ b/tests/test_synrm_5f_40s.py @@ -91,17 +91,18 @@ def test_Rstat_Lew_values(computation2): def test_extract_results_converts_units(computation2): computation2.set_variables(SimpleNamespace(variable_manager={}), Id1=1.0, Iq1=2.0, Id3=3.0, Iq3=4.0) - def data_real(expr): + def data_real(expr, convert_to_SI=False): + assert convert_to_SI, f"expected convert_to_SI=True for {expr}" letters = ["ABCDE".index(part[0]) + 1 for part in expr.split("Phase")[1:]] value = float(np.prod(letters) if letters else 7.0) - return np.full(2, value * 1e9 if expr.startswith("L(") else value) + return np.full(2, value) solutions = SimpleNamespace(data_real=data_real, primary_sweep_values=[0.0, 1.0]) out = computation2.extract_results(solutions) - theta_el = electrical_angle(np.full(2, 7.0), computation2.InitPos, computation2.geometry.PolePairs, computation2.RotSign, degrees=True) + theta_el = electrical_angle(np.full(2, 7.0), np.deg2rad(computation2.InitPos), computation2.geometry.PolePairs, computation2.RotSign, degrees=False) current_phases = np.stack([np.full(2, i + 1.0) for i in range(5)], axis=-1) expected_I_d1, _ = to_dq(current_phases, theta_el, harmonic=1) - assert out["I_d1"] == pytest.approx(expected_I_d1 / 1e3) + assert out["I_d1"] == pytest.approx(expected_I_d1) assert out["Moving1.Torque"] == pytest.approx(np.full(2, 7.0))