diff --git a/bluemira/magnets/cable.py b/bluemira/magnets/cable.py index 98df8d8113..78e2d5c1cd 100644 --- a/bluemira/magnets/cable.py +++ b/bluemira/magnets/cable.py @@ -6,52 +6,34 @@ """Cable class""" +from __future__ import annotations + from abc import ABC, abstractmethod -from collections.abc import Callable -from typing import Any +from typing import TYPE_CHECKING, Any import matplotlib.pyplot as plt import numpy as np from matproplib import OperationalConditions from scipy.integrate import solve_ivp -from scipy.optimize import minimize_scalar - -from bluemira.base.look_and_feel import bluemira_error, bluemira_print, bluemira_warn -from bluemira.magnets.registry import RegistrableMeta -from bluemira.magnets.strand import ( - Strand, - SuperconductingStrand, - create_strand_from_dict, -) -from bluemira.magnets.utils import parall_r, serie_r -# ------------------------------------------------------------------------------ -# Global Registries -# ------------------------------------------------------------------------------ -CABLE_REGISTRY = {} +from bluemira.base.look_and_feel import bluemira_debug +from bluemira.magnets.utils import reciprocal_summation, summation +if TYPE_CHECKING: + from collections.abc import Callable -# ------------------------------------------------------------------------------ -# Cable Class -# ------------------------------------------------------------------------------ + from bluemira.magnets.strand import Strand, SuperconductingStrand -class ABCCable(ABC, metaclass=RegistrableMeta): +class ABCCable(ABC): """ Abstract base class for superconducting cables. Defines the general structure and common methods for cables - composed of superconducting and stabilizer strands. + composed of superconducting and stabiliser strands. - Notes - ----- - - This class is abstract and cannot be instantiated directly. - - Subclasses must define `dx`, `dy`, `Kx`, `Ky`, and `from_dict`. """ - _registry_ = CABLE_REGISTRY - _name_in_registry_: str | None = None # Abstract base classes should NOT register - def __init__( self, sc_strand: SuperconductingStrand, @@ -59,9 +41,10 @@ def __init__( n_sc_strand: int, n_stab_strand: int, d_cooling_channel: float, - void_fraction: float = 0.725, - cos_theta: float = 0.97, + void_fraction: float, + cos_theta: float, name: str = "Cable", + **props, ): """ Representation of a cable. Only the x-dimension of the cable is given as @@ -73,42 +56,64 @@ def __init__( Parameters ---------- - sc_strand : SuperconductingStrand + sc_strand: The superconducting strand. - stab_strand : Strand - The stabilizer strand. - n_sc_strand : int + stab_strand: + The stabiliser strand. + n_sc_strand: Number of superconducting strands. - n_stab_strand : int - Number of stabilizing strands. - d_cooling_channel : float + n_stab_strand: + Number of stabilising strands. + d_cooling_channel: Diameter of the cooling channel [m]. - void_fraction : float + void_fraction: Ratio of material volume to total volume [unitless]. - cos_theta : float + cos_theta: Correction factor for twist in the cable layout. - name : str + name: Identifier for the cable instance. - """ - # initialize private variables - self._d_cooling_channel = None - self._void_fraction = None - self._n_sc_strand = None - self._n_stab_strand = None - self._cos_theta = None - self._shape = None + Raises + ------ + ValueError + If E not defined on the class and not passed in as a kwarg + """ # assign # Setting self.name triggers automatic instance registration self.name = name self.sc_strand = sc_strand self.stab_strand = stab_strand - self.void_fraction = void_fraction - self.d_cooling_channel = d_cooling_channel self.n_sc_strand = n_sc_strand self.n_stab_strand = n_stab_strand + self.d_cooling_channel = d_cooling_channel + self.void_fraction = void_fraction self.cos_theta = cos_theta + youngs_modulus: Callable[[Any, OperationalConditions], float] | float | None = ( + props.pop("E", None) + ) + + def ym(op_cond): + raise NotImplementedError("E for Cable is not implemented.") + + if "E" not in vars(type(self)): + if youngs_modulus is None: + youngs_modulus = ym + + self.E = ( + youngs_modulus + if callable(youngs_modulus) + else lambda op_cond, v=youngs_modulus: youngs_modulus # noqa: ARG005 + ) + elif youngs_modulus is not None: + bluemira_debug("E already defined in class, ignoring") + + for k, v in props.items(): + setattr(self, k, v if callable(v) else lambda *arg, v=v, **kwargs: v) # noqa: ARG005 + self._props = list(props.keys()) + ( + [] if "E" in vars(type(self)) or youngs_modulus == ym else ["E"] + ) + @property @abstractmethod def dx(self): @@ -126,105 +131,13 @@ def aspect_ratio(self): """ return self.dx / self.dy - @property - def n_sc_strand(self): - """Number of superconducting strands""" - return self._n_sc_strand - - @n_sc_strand.setter - def n_sc_strand(self, value: int): - """ - Set the number of superconducting strands. - - Raises - ------ - ValueError - If the value is not positive. - """ - if value <= 0: - msg = f"The number of superconducting strands must be positive, got {value}" - bluemira_error(msg) - raise ValueError(msg) - self._n_sc_strand = int(np.ceil(value)) - - @property - def n_stab_strand(self): - """Number of stabilizing strands""" - return self._n_stab_strand - - @n_stab_strand.setter - def n_stab_strand(self, value: int): - """ - Set the number of stabilizer strands. - - Raises - ------ - ValueError - If the value is negative. - """ - if value < 0: - msg = f"The number of stabilizing strands must be positive, got {value}" - bluemira_error(msg) - raise ValueError(msg) - self._n_stab_strand = int(np.ceil(value)) - - @property - def d_cooling_channel(self): - """Diameter of the cooling channel [m].""" - return self._d_cooling_channel - - @d_cooling_channel.setter - def d_cooling_channel(self, value: float): - """ - Set the cooling channel diameter. - - Raises - ------ - ValueError - If the value is negative. - """ - if value < 0: - msg = f"diameter of the cooling channel must be positive, got {value}" - bluemira_error(msg) - raise ValueError(msg) - - self._d_cooling_channel = value - - @property - def void_fraction(self): - """Void fraction of the cable.""" - return self._void_fraction - - @void_fraction.setter - def void_fraction(self, value: float): - if value < 0 or value > 1: - msg = f"void_fraction must be between 0 and 1, got {value}" - bluemira_error(msg) - raise ValueError(msg) - - self._void_fraction = value - - @property - def cos_theta(self): - """Correction factor for strand orientation (twist).""" - return self._cos_theta - - @cos_theta.setter - def cos_theta(self, value: float): - if value <= 0 or value > 1: - msg = f"cos theta must be in the interval ]0, 1], got {value}" - bluemira_error(msg) - raise ValueError(msg) - - self._cos_theta = value - def rho(self, op_cond: OperationalConditions): """ Compute the average mass density of the cable [kg/m³]. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. @@ -234,30 +147,31 @@ def rho(self, op_cond: OperationalConditions): Averaged mass density in kg/m³. """ return ( - self.sc_strand.rho(op_cond) * self.area_sc - + self.stab_strand.rho(op_cond) * self.area_stab - ) / (self.area_sc + self.area_stab) + self.sc_strand.rho(op_cond) * self.area_sc_region + + self.stab_strand.rho(op_cond) * self.area_stab_region + ) / (self.area_sc_region + self.area_stab_region) - def erho(self, op_cond: OperationalConditions): + def erho(self, op_cond: OperationalConditions) -> float: """ Computes the cable's equivalent resistivity considering the resistance of its strands in parallel. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float [Ohm m] + : + resistivity [Ohm m] """ resistances = np.array([ - self.sc_strand.erho(op_cond) / self.area_sc, - self.stab_strand.erho(op_cond) / self.area_stab, + self.sc_strand.erho(op_cond) / self.area_sc_region, + self.stab_strand.erho(op_cond) / self.area_stab_region, ]) - res_tot = parall_r(resistances) + res_tot = reciprocal_summation(resistances) return res_tot * self.area def Cp(self, op_cond: OperationalConditions): # noqa: N802 @@ -267,66 +181,49 @@ def Cp(self, op_cond: OperationalConditions): # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float [J/K/m] + : + Specific heat capacity [J/K/m] """ weighted_specific_heat = np.array([ - self.sc_strand.Cp(op_cond) * self.area_sc * self.sc_strand.rho(op_cond), + self.sc_strand.Cp(op_cond) + * self.area_sc_region + * self.sc_strand.rho(op_cond), self.stab_strand.Cp(op_cond) - * self.area_stab + * self.area_stab_region * self.stab_strand.rho(op_cond), ]) - return serie_r(weighted_specific_heat) / ( - self.area_sc * self.sc_strand.rho(op_cond) - + self.area_stab * self.stab_strand.rho(op_cond) + return summation(weighted_specific_heat) / ( + self.area_sc_region * self.sc_strand.rho(op_cond) + + self.area_stab_region * self.stab_strand.rho(op_cond) ) @property - def area_stab(self): - """Area of the stabilizer region""" + def area_stab_region(self) -> float: + """Area of the stabiliser region""" return self.stab_strand.area * self.n_stab_strand @property - def area_sc(self): + def area_sc_region(self) -> float: """Area of the superconductor region""" return self.sc_strand.area * self.n_sc_strand @property - def area_cc(self): + def area_cooling_channel(self) -> float: """Area of the cooling channel""" return self.d_cooling_channel**2 / 4 * np.pi @property - def area(self): + def area(self) -> float: """Area of the cable considering the void fraction""" return ( - self.area_sc + self.area_stab - ) / self.void_fraction / self.cos_theta + self.area_cc - - def E(self, op_cond: OperationalConditions): # noqa: N802 - """ - Return the effective Young's modulus of the cable [Pa]. - - This is a default placeholder implementation in the base class. - Subclasses may use `kwargs` to modify behavior. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Default Young's modulus (0). - """ - raise NotImplementedError("E for Cable is not implemented.") + self.area_sc_region + self.area_stab_region + ) / self.void_fraction / self.cos_theta + self.area_cooling_channel def _heat_balance_model_cable( self, @@ -334,29 +231,29 @@ def _heat_balance_model_cable( temperature: float, B_fun: Callable, I_fun: Callable, # noqa: N803 - ): + ) -> float: """ Calculate the derivative of temperature (dT/dt) for a 0D heat balance problem. Parameters ---------- - t : float - The current time in seconds. - temperature : float - The current temperature in Celsius. - B_fun : Callable - The magnetic field [T] as time function - I_fun : Callable - The current [A] flowing through the conductor as time function + t: + The current time in seconds. + temperature: + The current temperature in Celsius. + B_fun: + The magnetic field [T] as time function + I_fun: + The current [A] flowing through the conductor as time function Returns ------- - dTdt : float - The derivative of temperature with respect to time (dT/dt). + : + The derivative of temperature with respect to time (dT/dt). """ # Calculate the rate of heat generation (Joule dissipation) if isinstance(temperature, np.ndarray): - temperature = temperature[0] + temperature = temperature.item() op_cond = OperationalConditions(temperature=temperature, magnetic_field=B_fun(t)) @@ -375,8 +272,8 @@ def _temperature_evolution( t0: float, tf: float, initial_temperature: float, - B_fun: Callable, - I_fun: Callable, # noqa: N803 + B_fun: Callable[[float], float], + I_fun: Callable[[float], float], # noqa: N803 ): solution = solve_ivp( self._heat_balance_model_cable, @@ -391,225 +288,68 @@ def _temperature_evolution( return solution - def optimize_n_stab_ths( + def final_temperature_difference( self, + n_stab: int, t0: float, tf: float, initial_temperature: float, target_temperature: float, - B_fun: Callable, - I_fun: Callable, # noqa: N803 - bounds: np.ndarray = None, - *, - show: bool = False, - ): + B_fun: Callable[[float], float], + I_fun: Callable[[float], float], # noqa: N803 + ) -> float: """ - Optimize the number of stabilizer strand in the superconducting cable using a - 0-D hot spot criteria. + Compute the absolute temperature difference at final time between the + simulated and target temperatures. + + This method modifies the private attribute `_n_stab_strand` to update the + cable configuration, simulates the temperature evolution over time, and + returns the absolute difference between the final temperature and the + specified target. Parameters ---------- + n_stab: + Number of stabiliser strands to set temporarily for this simulation. t0: - Initial time [s]. + Initial time of the simulation [s]. tf: - Final time [s]. + Final time of the simulation [s]. initial_temperature: - Temperature [K] at initial time. + Temperature at the start of the simulation [K]. target_temperature: - Target temperature [K] at final time. - B_fun : - Magnetic field [T] as a time-dependent function. - I_fun : - Current [A] as a time-dependent function. - bounds: - Lower and upper limits for the number of stabilizer strands. - show: - If True, the behavior of temperature over time is plotted. + Desired temperature at the end of the simulation [K]. + B_fun: + Magnetic field as a time-dependent function [T]. + I_fun: + Current as a time-dependent function [A]. Returns ------- - result : scipy.optimize.OptimizeResult - The result of the optimization process. - - Raises - ------ - ValueError - If the optimization process does not converge. + : + Absolute difference between the simulated final temperature and the + target temperature [K]. Notes ----- - - The number of stabilizer strands in the cable is modified directly. - - Cooling material contribution is neglected when applying the hot spot criteria. - """ - - def final_temperature_difference( - n_stab: int, - t0: float, - tf: float, - initial_temperature: float, - target_temperature: float, - B_fun: Callable, - I_fun: Callable, # noqa: N803 - ): - """ - Compute the absolute temperature difference at final time between the - simulated and target temperatures. - - This method modifies the private attribute `_n_stab_strand` to update the - cable configuration, simulates the temperature evolution over time, and - returns the absolute difference between the final temperature and the - specified target. - - Parameters - ---------- - n_stab : int - Number of stabilizer strands to set temporarily for this simulation. - t0 : float - Initial time of the simulation [s]. - tf : float - Final time of the simulation [s]. - initial_temperature : float - Temperature at the start of the simulation [K]. - target_temperature : float - Desired temperature at the end of the simulation [K]. - B_fun : Callable - Magnetic field as a time-dependent function [T]. - I_fun : Callable - Current as a time-dependent function [A]. - - Returns - ------- - float - Absolute difference between the simulated final temperature and the - target temperature [K]. - - Notes - ----- - - This method is typically used as a cost function for optimization routines - (e.g., minimizing the temperature error by tuning `n_stab`). - - It modifies the internal state `self._n_stab_strand`, which may affect - subsequent evaluations unless restored. - """ - self._n_stab_strand = n_stab - - solution = self._temperature_evolution( - t0=t0, - tf=tf, - initial_temperature=initial_temperature, - B_fun=B_fun, - I_fun=I_fun, - ) - final_temperature = float(solution.y[0][-1]) - # diff = abs(final_temperature - target_temperature) - return abs(final_temperature - target_temperature) - - method = None - if bounds is not None: - method = "bounded" - - result = minimize_scalar( - fun=final_temperature_difference, - args=(t0, tf, initial_temperature, target_temperature, B_fun, I_fun), - bounds=bounds, - method=method, + - This method is typically used as a cost function for optimisation routines + (e.g., minimising the temperature error by tuning `n_stab`). + - It modifies the internal state `self._n_stab_strand`, which may affect + subsequent evaluations unless restored. + """ + self.n_stab_strand = n_stab + + solution = self._temperature_evolution( + t0=t0, + tf=tf, + initial_temperature=initial_temperature, + B_fun=B_fun, + I_fun=I_fun, ) + final_temperature = float(solution.y[0][-1]) + # diff = abs(final_temperature - target_temperature) + return abs(final_temperature - target_temperature) - if not result.success: - raise ValueError( - "n_stab optimization did not converge. Check your input parameters " - "or initial bracket." - ) - - # Here we re-ensure the n_stab_strand to be an integer - self.n_stab_strand = self._n_stab_strand - - solution = self._temperature_evolution(t0, tf, initial_temperature, B_fun, I_fun) - final_temperature = solution.y[0][-1] - - if final_temperature > target_temperature: - bluemira_error( - f"Final temperature ({final_temperature:.2f} K) exceeds target " - f"temperature " - f"({target_temperature} K) even with maximum n_stab = " - f"{self.n_stab_strand}." - ) - raise ValueError( - "Optimization failed to keep final temperature ≤ target. " - "Try increasing the upper bound of n_stab or adjusting cable parameters." - ) - bluemira_print(f"Optimal n_stab: {self.n_stab_strand}") - bluemira_print( - f"Final temperature with optimal n_stab: {final_temperature:.2f} Kelvin" - ) - - if show: - _, (ax_temp, ax_ib) = plt.subplots(2, 1, figsize=(8, 8), sharex=True) - - # --- Plot Temperature Evolution --- - ax_temp.plot(solution.t, solution.y[0], "r*", label="Simulation points") - time_steps = np.linspace(t0, tf, 100) - ax_temp.plot( - time_steps, solution.sol(time_steps)[0], "b", label="Interpolated curve" - ) - ax_temp.grid(visible=True) - ax_temp.set_ylabel("Temperature [K]", fontsize=10) - ax_temp.set_title("Quench temperature evolution", fontsize=11) - ax_temp.legend(fontsize=9) - - ax_temp.tick_params(axis="y", labelcolor="k", labelsize=9) - - # Insert text box with additional info - info_text = ( - f"Target T: {target_temperature:.2f} K\n" - f"Initial T: {initial_temperature:.2f} K\n" - f"SC Strand: {self.sc_strand.name}\n" - f"n. sc. strand = {self.n_sc_strand}\n" - f"Stab. strand = {self.stab_strand.name}\n" - f"n. stab. strand = {self.n_stab_strand}\n" - ) - props = {"boxstyle": "round", "facecolor": "white", "alpha": 0.8} - ax_temp.text( - 0.65, - 0.5, - info_text, - transform=ax_temp.transAxes, - fontsize=9, - verticalalignment="top", - bbox=props, - ) - - # --- Plot I_fun(t) and B_fun(t) --- - time_steps_fine = np.linspace(t0, tf, 300) - I_values = [I_fun(t) for t in time_steps_fine] # noqa: N806 - B_values = [B_fun(t) for t in time_steps_fine] - - ax_ib.plot(time_steps_fine, I_values, "g", label="Current [A]") - ax_ib.set_ylabel("Current [A]", color="g", fontsize=10) - ax_ib.tick_params(axis="y", labelcolor="g", labelsize=9) - ax_ib.grid(visible=True) - - ax_ib_right = ax_ib.twinx() - ax_ib_right.plot( - time_steps_fine, B_values, "m--", label="Magnetic field [T]" - ) - ax_ib_right.set_ylabel("Magnetic field [T]", color="m", fontsize=10) - ax_ib_right.tick_params(axis="y", labelcolor="m", labelsize=9) - - # Labels - ax_ib.set_xlabel("Time [s]", fontsize=10) - ax_ib.tick_params(axis="x", labelsize=9) - - # Combined legend for both sides - lines, labels = ax_ib.get_legend_handles_labels() - lines2, labels2 = ax_ib_right.get_legend_handles_labels() - ax_ib.legend(lines + lines2, labels + labels2, loc="best", fontsize=9) - - plt.tight_layout() - plt.show() - - return result - - # OD homogenized structural properties @abstractmethod def Kx(self, op_cond: OperationalConditions): # noqa: N802 """Total equivalent stiffness along x-axis""" @@ -618,32 +358,39 @@ def Kx(self, op_cond: OperationalConditions): # noqa: N802 def Ky(self, op_cond: OperationalConditions): # noqa: N802 """Total equivalent stiffness along y-axis""" - def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): + def plot( + self, + xc: float = 0, + yc: float = 0, + *, + show: bool = False, + ax: plt.Axes | None = None, + ): """ Plot a schematic view of the cable cross-section. - This method visualizes the outer shape of the cable and the cooling channel, + This method visualises the outer shape of the cable and the cooling channel, assuming a rectangular or elliptical layout based on `dx`, `dy`, and `d_cooling_channel`. It draws the cable centered at (xc, yc) within the current coordinate system. Parameters ---------- - xc : float, optional + xc: x-coordinate of the cable center in the plot [m]. Default is 0. - yc : float, optional + yc: y-coordinate of the cable center in the plot [m]. Default is 0. - show : bool, optional + show: If True, the plot is rendered immediately with `plt.show()`. Default is False. - ax : matplotlib.axes.Axes or None, optional + ax: The matplotlib Axes object to draw on. If None, a new figure and Axes are created internally. Returns ------- - matplotlib.axes.Axes - The Axes object with the cable plot, which can be further customized + : + The Axes object with the cable plot, which can be further customised or saved. Notes @@ -658,6 +405,7 @@ def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): _, ax = plt.subplots() pc = np.array([xc, yc]) + a = self.dx / 2 b = self.dy / 2 @@ -683,16 +431,16 @@ def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): plt.show() return ax - def __str__(self): + def __str__(self) -> str: """ Return a human-readable summary of the cable configuration. Includes geometric properties, void and twist factors, and a string - representation of both the superconducting and stabilizer strands. + representation of both the superconducting and stabiliser strands. Returns ------- - str + : A formatted multiline string describing the cable. """ return ( @@ -712,9 +460,9 @@ def __str__(self): f"n stab strand: {self.n_stab_strand}" ) - def to_dict(self) -> dict: + def to_dict(self, op_cond) -> dict[str, str | float | int | dict[str, Any]]: """ - Serialize the cable instance to a dictionary. + Serialise the cable instance to a dictionary. Returns ------- @@ -722,9 +470,6 @@ def to_dict(self) -> dict: Dictionary containing cable and strand configuration. """ return { - "name_in_registry": getattr( - self, "_name_in_registry_", self.__class__.__name__ - ), "name": self.name, "n_sc_strand": self.n_sc_strand, "n_stab_strand": self.n_stab_strand, @@ -733,70 +478,9 @@ def to_dict(self) -> dict: "cos_theta": self.cos_theta, "sc_strand": self.sc_strand.to_dict(), "stab_strand": self.stab_strand.to_dict(), + **{k: getattr(self, k)(op_cond) for k in self._props}, } - @classmethod - def from_dict( - cls, - cable_dict: dict[str, Any], - name: str | None = None, - ) -> "ABCCable": - """ - Deserialize a cable instance from a dictionary. - - Parameters - ---------- - cls : type - Class to instantiate (Cable or subclass). - cable_dict : dict - Dictionary containing serialized cable data. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - ABCCable - Instantiated cable object. - - Raises - ------ - ValueError - If name_in_registry mismatch or duplicate instance name. - """ - name_in_registry = cable_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - # Deserialize strands - sc_strand_data = cable_dict["sc_strand"] - if isinstance(sc_strand_data, Strand): - sc_strand = sc_strand_data - else: - sc_strand = create_strand_from_dict(strand_dict=sc_strand_data) - - stab_strand_data = cable_dict["stab_strand"] - if isinstance(stab_strand_data, Strand): - stab_strand = stab_strand_data - else: - stab_strand = create_strand_from_dict(strand_dict=stab_strand_data) - - return cls( - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=cable_dict["n_sc_strand"], - n_stab_strand=cable_dict["n_stab_strand"], - d_cooling_channel=cable_dict["d_cooling_channel"], - void_fraction=cable_dict.get("void_fraction", 0.725), - cos_theta=cable_dict.get("cos_theta", 0.97), - name=name or cable_dict.get("name"), - ) - class RectangularCable(ABCCable): """ @@ -806,19 +490,18 @@ class RectangularCable(ABCCable): the total area and x-dimension. """ - _name_in_registry_ = "RectangularCable" - def __init__( self, - dx: float, sc_strand: SuperconductingStrand, stab_strand: Strand, n_sc_strand: int, n_stab_strand: int, d_cooling_channel: float, - void_fraction: float = 0.725, - cos_theta: float = 0.97, + void_fraction: float, + cos_theta: float, + dx: float, name: str = "RectangularCable", + **props, ): """ Representation of a cable. Only the x-dimension of the cable is given as @@ -830,24 +513,26 @@ def __init__( Parameters ---------- - dx : float - Cable width in the x-direction [m]. - sc_strand : SuperconductingStrand + sc_strand: Superconducting strand. - stab_strand : Strand - Stabilizer strand. - n_sc_strand : int + stab_strand: + Stabiliser strand. + n_sc_strand: Number of superconducting strands. - n_stab_strand : int - Number of stabilizer strands. - d_cooling_channel : float - Cooling channel diameter [m]. - void_fraction : float, optional - Void fraction (material_volume / total_volume). - cos_theta : float, optional - Correction factor for strand twist. - name : str, optional - Name of the cable. + n_stab_strand: + Number of stabilising strands. + d_cooling_channel: + Diameter of the cooling channel [m]. + void_fraction: + Ratio of material volume to total volume [unitless]. + cos_theta: + Correction factor for twist in the cable layout. + dx: + Cable half-width in the x-direction [m]. + name: + Name of the cable + props: + extra properties """ super().__init__( sc_strand=sc_strand, @@ -858,278 +543,83 @@ def __init__( void_fraction=void_fraction, cos_theta=cos_theta, name=name, + **props, ) - - # initialize private variables - self._dx = None - - # assign - self.dx = dx + self._dx = dx @property - def dx(self): + def dx(self) -> float: """Cable dimension in the x direction [m]""" return self._dx - @dx.setter - def dx(self, value: float): - """ - Set cable width in x-direction. - - Raises - ------ - ValueError - If value is not positive. - """ - if value <= 0: - msg = "dx must be positive" - bluemira_error(msg) - raise ValueError(msg) - self._dx = value - @property - def dy(self): + def dy(self) -> float: """Cable dimension in the y direction [m]""" return self.area / self.dx - # Decide if this function shall be a setter. - # Defined as "normal" function to underline that it modifies dx. - def set_aspect_ratio(self, value: float) -> None: + @ABCCable.aspect_ratio.setter + def aspect_ratio(self, value: float): """Modify dx in order to get the given aspect ratio""" - self.dx = np.sqrt(value * self.area) + self._dx = np.sqrt(value * self.area) - # OD homogenized structural properties - def Kx(self, op_cond: OperationalConditions): # noqa: N802 + def Kx(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the total equivalent stiffness along the x-axis. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float - Homogenized stiffness in the x-direction [Pa]. + : + Homogenised stiffness in the x-direction [Pa]. """ - return self.E(op_cond) * self.dy / self.dx + return self.E(op_cond) / self.aspect_ratio - def Ky(self, op_cond: OperationalConditions): # noqa: N802 + def Ky(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the total equivalent stiffness along the y-axis. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float - Homogenized stiffness in the y-direction [Pa]. + : + Homogenised stiffness in the y-direction [Pa]. """ - return self.E(op_cond) * self.dx / self.dy + return self.E(op_cond) * self.aspect_ratio - def to_dict(self) -> dict: + def to_dict(self, op_cond) -> dict[str, Any]: """ - Serialize the rectangular cable into a dictionary. + Serialise the rectangular cable into a dictionary. Returns ------- - dict + : Dictionary including rectangular cable parameters. """ - data = super().to_dict() + data = super().to_dict(op_cond) data.update({ "dx": self.dx, "aspect_ratio": self.aspect_ratio, }) return data - @classmethod - def from_dict( - cls, - cable_dict: dict[str, Any], - name: str | None = None, - ) -> "RectangularCable": - """ - Deserialize a RectangularCable from a dictionary. - - Behavior: - - If both 'dx' and 'aspect_ratio' are provided, a warning is issued and - aspect_ratio is applied. - - If only 'aspect_ratio' is provided, dx and dy are calculated accordingly. - - If only 'dx' is provided, it is used as-is. - - If neither is provided, raises a ValueError. - - Parameters - ---------- - cls : type - Class to instantiate (Cable or subclass). - cable_dict : dict - Dictionary containing serialized cable data. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - RectangularCable - Instantiated rectangular cable object. - - Raises - ------ - ValueError - If neither 'dx' nor 'aspect_ratio' is provided. - """ - name_in_registry = cable_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - # Deserialize strands - sc_strand_data = cable_dict["sc_strand"] - if isinstance(sc_strand_data, Strand): - sc_strand = sc_strand_data - else: - sc_strand = create_strand_from_dict(strand_dict=sc_strand_data) - - stab_strand_data = cable_dict["stab_strand"] - if isinstance(stab_strand_data, Strand): - stab_strand = stab_strand_data - else: - stab_strand = create_strand_from_dict(strand_dict=stab_strand_data) - - # Geometry parameters - dx = cable_dict.get("dx") - aspect_ratio = cable_dict.get("aspect_ratio") - - if dx is not None and aspect_ratio is not None: - bluemira_warn( - "Both 'dx' and 'aspect_ratio' specified. Aspect ratio will override dx " - "after creation." - ) - - if aspect_ratio is not None and dx is None: - # Default dx if only aspect ratio is provided. It will be recalculated at - # the end when set_aspect_ratio is called - dx = 0.01 - if dx is None: - raise ValueError( - "Serialized RectangularCable must include at least 'dx' or " - "'aspect_ratio'." - ) - - # Base cable parameters - n_sc_strand = cable_dict["n_sc_strand"] - n_stab_strand = cable_dict["n_stab_strand"] - d_cooling_channel = cable_dict["d_cooling_channel"] - void_fraction = cable_dict.get("void_fraction", 0.725) - cos_theta = cable_dict.get("cos_theta", 0.97) - - # Create cable - cable = cls( - dx=dx, - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=n_sc_strand, - n_stab_strand=n_stab_strand, - d_cooling_channel=d_cooling_channel, - void_fraction=void_fraction, - cos_theta=cos_theta, - name=name or cable_dict.get("name"), - ) - - # Adjust aspect ratio if needed - if aspect_ratio is not None: - cable.set_aspect_ratio(aspect_ratio) - - return cable - - -class DummyRectangularCableHTS(RectangularCable): - """ - Dummy rectangular cable with young's moduli set to 120 GPa. - """ - - _name_in_registry_ = "DummyRectangularCableHTS" - - def __init__(self, *args, **kwargs): - kwargs.setdefault("name", "DummyRectangularCableHTS") - super().__init__(*args, **kwargs) - - def E(self, op_cond: OperationalConditions): # noqa: N802, PLR6301, ARG002 - """ - Return the Young's modulus of the cable material. - - This is a constant value specific to the implementation. Subclasses may override - this method to provide a temperature- or field-dependent modulus. The `kwargs` - parameter is unused here but retained for interface consistency. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Young's modulus in Pascals [Pa]. - """ - return 120e9 - - -class DummyRectangularCableLTS(RectangularCable): - """ - Dummy square cable with young's moduli set to 0.1 GPa - """ - - _name_in_registry_ = "DummyRectangularCableLTS" - - def __init__(self, *args, **kwargs): - kwargs.setdefault("name", "DummyRectangularCableLTS") - super().__init__(*args, **kwargs) - - def E(self, op_cond): # noqa: N802, PLR6301, ARG002 - """ - Return the Young's modulus of the cable material. - - This implementation returns a fixed value (0.1 GPa). Subclasses may override - this method with more sophisticated behavior. `kwargs` are included for - compatibility but not used in this implementation. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Young's modulus in Pascals [Pa]. - """ - return 0.1e9 - - -class SquareCable(ABCCable): +class SquareCable(RectangularCable): """ Cable with a square cross-section. Both dx and dy are derived from the total cross-sectional area. """ - _name_in_registry_ = "SquareCable" - def __init__( self, sc_strand: SuperconductingStrand, @@ -1137,9 +627,10 @@ def __init__( n_sc_strand: int, n_stab_strand: int, d_cooling_channel: float, - void_fraction: float = 0.725, - cos_theta: float = 0.97, + void_fraction: float, + cos_theta: float, name: str = "SquareCable", + **props, ): """ Representation of a square cable. @@ -1154,17 +645,17 @@ def __init__( sc_strand: strand of the superconductor stab_strand: - strand of the stabilizer - d_cooling_channel: - diameter of the cooling channel + strand of the stabiliser n_sc_strand: - number of superconducting strands + Number of superconducting strands. n_stab_strand: - number of stabilizer strands + Number of stabilising strands. + d_cooling_channel: + Diameter of the cooling channel [m]. void_fraction: - void fraction defined as material_volume/total_volume + Ratio of material volume to total volume [unitless]. cos_theta: - corrective factor that consider the twist of the cable + Correction factor for twist in the cable layout. name: cable string identifier @@ -1181,186 +672,39 @@ def __init__( void_fraction=void_fraction, cos_theta=cos_theta, name=name, + **props, ) @property - def dx(self): + def dx(self) -> float: """Cable dimension in the x direction [m]""" return np.sqrt(self.area) @property - def dy(self): + def dy(self) -> float: """Cable dimension in the y direction [m]""" return self.dx - # OD homogenized structural properties - def Kx(self, op_cond: OperationalConditions): # noqa: N802 - """ - Compute the total equivalent stiffness along the x-axis. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Homogenized stiffness in the x-direction [Pa]. - """ - return self.E(op_cond) * self.dy / self.dx - - def Ky(self, op_cond: OperationalConditions): # noqa: N802 - """ - Compute the total equivalent stiffness along the y-axis. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Homogenized stiffness in the y-direction [Pa]. - """ - return self.E(op_cond) * self.dx / self.dy - - def to_dict(self) -> dict: - """ - Serialize the SquareCable. - - Returns - ------- - dict - Serialized dictionary. - """ - return super().to_dict() - - @classmethod - def from_dict( - cls, - cable_dict: dict[str, Any], - name: str | None = None, - ) -> "SquareCable": - """ - Deserialize a SquareCable from a dictionary. - - Parameters - ---------- - cls : type - Class to instantiate (Cable or subclass). - cable_dict : dict - Dictionary containing serialized cable data. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - SquareCable - Instantiated square cable. - - Raises - ------ - ValueError - If unique_name is False and a duplicate name is detected in the instance - cache. - """ - name_in_registry = cable_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - sc_strand = create_strand_from_dict(strand_dict=cable_dict["sc_strand"]) - stab_strand = create_strand_from_dict(strand_dict=cable_dict["stab_strand"]) - - return cls( - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=cable_dict["n_sc_strand"], - n_stab_strand=cable_dict["n_stab_strand"], - d_cooling_channel=cable_dict["d_cooling_channel"], - void_fraction=cable_dict.get("void_fraction", 0.725), - cos_theta=cable_dict.get("cos_theta", 0.97), - name=name or cable_dict.get("name"), - ) - - -class DummySquareCableHTS(SquareCable): - """ - Dummy square cable with Young's modulus set to 120 GPa. - """ - - _name_in_registry_ = "DummySquareCableHTS" - - def __init__(self, *args, **kwargs): - kwargs.setdefault("name", "DummySquareCableHTS") - super().__init__(*args, **kwargs) - - def E(self, op_cond: OperationalConditions): # noqa: N802, PLR6301, ARG002 - """ - Return the Young's modulus for the HTS dummy cable. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Young's modulus in Pascals [Pa]. + @property + def aspect_ratio(self): """ - return 120e9 - - -class DummySquareCableLTS(SquareCable): - """ - Dummy square cable with Young's modulus set to 0.1 GPa. - """ - - _name_in_registry_ = "DummySquareCableLTS" - - def __init__(self, *args, **kwargs): - kwargs.setdefault("name", "DummySquareCableLTS") - super().__init__(*args, **kwargs) - - def E(self, op_cond: OperationalConditions): # noqa: N802, PLR6301, ARG002 + Compute the aspect ratio of the cable cross-section. """ - Return the Young's modulus for the LTS dummy cable. + return 1 - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Young's modulus in Pascals [Pa]. - """ - return 0.1e9 + @aspect_ratio.setter + def aspect_ratio(self, _: Any): + raise AttributeError(f"Aspect Ratio cannot be set on {type(self)}") class RoundCable(ABCCable): """ A cable with round cross-section for superconducting applications. - This cable type includes superconducting and stabilizer strands arranged + This cable type includes superconducting and stabiliser strands arranged around a central cooling channel. """ - _name_in_registry_ = "RoundCable" - def __init__( self, sc_strand: SuperconductingStrand, @@ -1368,9 +712,10 @@ def __init__( n_sc_strand: int, n_stab_strand: int, d_cooling_channel: float, - void_fraction: float = 0.725, - cos_theta: float = 0.97, + void_fraction: float, + cos_theta: float, name: str = "RoundCable", + **props, ): """ Representation of a round cable @@ -1380,17 +725,17 @@ def __init__( sc_strand: strand of the superconductor stab_strand: - strand of the stabilizer - d_cooling_channel: - diameter of the cooling channel + strand of the stabiliser n_sc_strand: - number of superconducting strands + Number of superconducting strands. n_stab_strand: - number of stabilizer strands + Number of stabilising strands. + d_cooling_channel: + Diameter of the cooling channel [m]. void_fraction: - void fraction defined as material_volume/total_volume + Ratio of material volume to total volume [unitless]. cos_theta: - corrective factor that consider the twist of the cable + Correction factor for twist in the cable layout. name: cable string identifier """ @@ -1403,43 +748,44 @@ def __init__( void_fraction=void_fraction, cos_theta=cos_theta, name=name, + **props, ) @property - def dx(self): + def dx(self) -> float: """Cable dimension in the x direction [m] (i.e. cable's diameter)""" return np.sqrt(self.area * 4 / np.pi) @property - def dy(self): + def dy(self) -> float: """Cable dimension in the y direction [m] (i.e. cable's diameter)""" return self.dx - # OD homogenized structural properties - # A structural analysis should be performed to check how much the rectangular - # approximation is fine also for the round cable. - def Kx(self, op_cond: OperationalConditions): # noqa: N802 + # TODO: A structural analysis should be performed to check how much the rectangular + # approximation is fine also for the round cable. + def Kx(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the equivalent stiffness of the cable along the x-axis. - This is a homogenized 1D structural property derived from the Young's modulus + This is a homogenised 1D structural property derived from the Young's modulus and the cable's geometry. The stiffness reflects the effective resistance to deformation in the x-direction. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Equivalent stiffness in the x-direction [Pa]. """ - return self.E(op_cond) * self.dy / self.dx + # TODO possible reason for floating point difference + return self.E(op_cond) - def Ky(self, op_cond: OperationalConditions): # noqa: N802 + def Ky(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the equivalent stiffness of the cable along the y-axis. @@ -1449,37 +795,45 @@ def Ky(self, op_cond: OperationalConditions): # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Equivalent stiffness in the y-direction [Pa]. """ - return self.E(op_cond) * self.dx / self.dy + # TODO possible reason for floating point difference + return self.E(op_cond) - def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): + def plot( + self, + xc: float = 0, + yc: float = 0, + *, + show: bool = False, + ax: plt.Axes | None = None, + ) -> plt.Axes: """ Schematic plot of the cable cross-section. Parameters ---------- - xc : float, optional + xc: x-coordinate of the cable center [m]. Default is 0. - yc : float, optional + yc: y-coordinate of the cable center [m]. Default is 0. - show : bool, optional + show: If True, the plot is displayed immediately using `plt.show()`. Default is False. - ax : matplotlib.axes.Axes or None, optional + ax: Axis to plot on. If None, a new figure and axis are created. Returns ------- - matplotlib.axes.Axes - The axis object containing the cable plot, useful for further customization + : + The axis object containing the cable plot, useful for further customisation or saving. """ if ax is None: @@ -1509,171 +863,3 @@ def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): if show: plt.show() return ax - - def to_dict(self) -> dict: - """ - Serialize the RoundCable. - - Returns - ------- - dict - Serialized dictionary. - """ - return super().to_dict() - - @classmethod - def from_dict( - cls, - cable_dict: dict[str, Any], - name: str | None = None, - ) -> "RoundCable": - """ - Deserialize a RoundCable from a dictionary. - - Parameters - ---------- - cls : type - Class to instantiate (Cable or subclass). - cable_dict : dict - Dictionary containing serialized cable data. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - RoundCable - Instantiated square cable. - - Raises - ------ - ValueError - If unique_name is False and a duplicate name is detected in the instance - cache. - """ - name_in_registry = cable_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - sc_strand = create_strand_from_dict(strand_dict=cable_dict["sc_strand"]) - stab_strand = create_strand_from_dict(strand_dict=cable_dict["stab_strand"]) - - return cls( - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=cable_dict["n_sc_strand"], - n_stab_strand=cable_dict["n_stab_strand"], - d_cooling_channel=cable_dict["d_cooling_channel"], - void_fraction=cable_dict.get("void_fraction", 0.725), - cos_theta=cable_dict.get("cos_theta", 0.97), - name=name or cable_dict.get("name"), - ) - - -class DummyRoundCableHTS(RoundCable): - """ - Dummy round cable with Young's modulus set to 120 GPa. - - This class provides a simplified round cable configuration for high-temperature - superconducting (HTS) analysis with a fixed stiffness value. - """ - - _name_in_registry_ = "DummyRoundCableHTS" - - def __init__(self, *args, **kwargs): - kwargs.setdefault("name", "DummyRoundCableHTS") - super().__init__(*args, **kwargs) - - def E(self, op_cond: OperationalConditions): # noqa: N802, PLR6301, ARG002 - """ - Return the Young's modulus for the HTS dummy round cable. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Young's modulus in Pascals [Pa]. - """ - return 120e9 - - -class DummyRoundCableLTS(RoundCable): - """ - Dummy round cable with Young's modulus set to 0.1 GPa. - - This class provides a simplified round cable configuration for low-temperature - superconducting (LTS) analysis with a fixed, softer stiffness value. - """ - - _name_in_registry_ = "DummyRoundCableLTS" - - def __init__(self, *args, **kwargs): - kwargs.setdefault("name", "DummyRoundCableLTS") - super().__init__(*args, **kwargs) - - def E(self, op_cond: OperationalConditions): # noqa: N802, PLR6301, ARG002 - """ - Return the Young's modulus for the LTS dummy round cable. - - Parameters - ---------- - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - - Returns - ------- - float - Young's modulus in Pascals [Pa]. - """ - return 0.1e9 - - -def create_cable_from_dict( - cable_dict: dict, - name: str | None = None, -): - """ - Factory function to create a Cable or its subclass from a serialized dictionary. - - Parameters - ---------- - cable_dict : dict - Dictionary with serialized cable data. Must include a 'name_in_registry' field. - name : str, optional - If given, overrides the name from the dictionary. - - Returns - ------- - ABCCable - Instantiated cable object. - - Raises - ------ - ValueError - If 'name_in_registry' is missing or no matching class is found. - """ - name_in_registry = cable_dict.get("name_in_registry") - if name_in_registry is None: - raise ValueError( - "Serialized cable dictionary must contain a 'name_in_registry' field." - ) - - cls = CABLE_REGISTRY.get(name_in_registry) - if cls is None: - raise ValueError( - f"No registered cable class with registration name '{name_in_registry}'. " - "Available classes are: " + ", ".join(CABLE_REGISTRY.keys()) - ) - - return cls.from_dict(name=name, cable_dict=cable_dict) diff --git a/bluemira/magnets/case_tf.py b/bluemira/magnets/case_tf.py index 7031801e8c..1f0610206d 100644 --- a/bluemira/magnets/case_tf.py +++ b/bluemira/magnets/case_tf.py @@ -7,7 +7,7 @@ """ Toroidal Field (TF) Coil 2D Case Class. -This class models and optimizes the cross-sectional layout of the inboard leg of a TF +This class models and optimises the cross-sectional layout of the inboard leg of a TF coil. It is designed to define and adjust the distribution of structural materials and winding pack arrangement to achieve optimal performance and mechanical robustness. @@ -16,235 +16,101 @@ - Focused on the two-dimensional analysis of the inboard leg. """ +from __future__ import annotations + import math from abc import ABC, abstractmethod +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any import matplotlib.pyplot as plt import numpy as np -from matproplib import OperationalConditions -from matproplib.material import Material -from scipy.optimize import minimize_scalar from bluemira.base.look_and_feel import ( bluemira_debug, bluemira_error, - bluemira_print, bluemira_warn, ) -from bluemira.magnets.registry import RegistrableMeta -from bluemira.magnets.utils import parall_k, serie_k -from bluemira.magnets.winding_pack import WindingPack, create_wp_from_dict +from bluemira.geometry.parameterisations import GeometryParameterisation +from bluemira.geometry.tools import make_polygon +from bluemira.magnets.utils import reciprocal_summation, summation +from bluemira.magnets.winding_pack import WindingPack +from bluemira.utilities.opt_variables import OptVariable, OptVariablesFrame, VarDictT, ov -# ------------------------------------------------------------------------------ -# Global Registries -# ------------------------------------------------------------------------------ -CASETF_REGISTRY = {} +if TYPE_CHECKING: + from matproplib import OperationalConditions + from matproplib.material import Material + from bluemira.geometry.wire import BluemiraWire -# ------------------------------------------------------------------------------ -# TFcoil cross section Geometry Base and Implementations -# ------------------------------------------------------------------------------ -class CaseGeometry(ABC): - """ - Abstract base class for TF case geometry profiles. - Provides access to radial dimensions and toroidal width calculations - as well as geometric plotting and area calculation interfaces. +def _dx_at_radius(radius: float, rad_theta: float) -> float: """ + Compute the toroidal width at a given radial position. - def __init__(self, Ri: float, Rk: float, theta_TF: float): # noqa: N803 - """ - Initialize the geometry base. - - Parameters - ---------- - Ri : float - External radius of the TF coil case [m]. - Rk : float - Internal radius of the TF coil case [m]. - theta_TF : float - Toroidal angular span of the TF coil [degrees]. - """ - self._Ri = None - self.Ri = Ri - - self._Rk = None - self.Rk = Rk - - self.theta_TF = theta_TF - - @property - def Ri(self) -> float: # noqa: N802 - """ - External (outermost) radius of the TF case at the top [m]. - - Returns - ------- - float - Outer radius measured from the machine center to the case outer wall [m]. - """ - return self._Ri - - @Ri.setter - def Ri(self, value: float): # noqa: N802 - """ - Set the external (outermost) radius of the TF case. - - Parameters - ---------- - value : float - Outer radius [m]. Must be a strictly positive number. - - Raises - ------ - ValueError - If the provided radius is not positive. - """ - if value <= 0: - raise ValueError("Ri must be positive.") - self._Ri = value - - @property - def Rk(self) -> float: # noqa: N802 - """ - Internal (innermost) radius of the TF case at the top [m]. - - Returns - ------- - float - Inner radius measured from the machine center to the case outer wall [m]. - """ - return self._Rk - - @Rk.setter - def Rk(self, value: float): # noqa: N802 - """ - Set the internal (innermost) radius of the TF case. - - Parameters - ---------- - value : float - Outer radius [m]. Must be a strictly positive number. - - Raises - ------ - ValueError - If the provided radius is not positive. - """ - if value < 0: - raise ValueError("Rk must be positive.") - self._Rk = value - - @property - def theta_TF(self) -> float: - """ - Toroidal angular span of the TF coil [degrees]. - - Returns - ------- - float - Toroidal angular span [°]. - """ - return self._theta_TF - - @theta_TF.setter - def theta_TF(self, value: float): - """ - Set the toroidal angular span and update the internal radian representation. - - Parameters - ---------- - value : float - New toroidal angular span [degrees]. - - Raises - ------ - ValueError - If the provided value is not within (0, 360] degrees. - """ - if not (0.0 < value <= 360.0): # noqa: PLR2004 - raise ValueError("theta_TF must be in the range (0, 360] degrees.") - self._theta_TF = value - self._rad_theta_TF = np.radians(value) - - @property - def rad_theta_TF(self): - """ + Parameters + ---------- + radius: + Radial position at which to compute the toroidal width [m]. + rad_theta: Toroidal angular span of the TF coil [radians]. - Returns - ------- - float - Toroidal aperture converted to radians. - """ - return self._rad_theta_TF - - def dx_at_radius(self, radius: float) -> float: - """ - Compute the toroidal width at a given radial position. - - Parameters - ---------- - radius : float - Radial position at which to compute the toroidal width [m]. - - Returns - ------- - float - Toroidal width [m] at the given radius. - """ - return 2 * radius * np.tan(self.rad_theta_TF / 2) - - @property - @abstractmethod - def area(self) -> float: - """ - Compute the cross-sectional area of the TF case. - - Returns - ------- - float - Cross-sectional area [m²] enclosed by the case geometry. - - Notes - ----- - Must be implemented by each specific geometry class. - """ - - @abstractmethod - def plot(self, ax=None, *, show: bool = False) -> plt.Axes: - """ - Plot the cross-sectional geometry of the TF case. + Returns + ------- + : + Toroidal width [m] at the given radius. + """ + return 2 * radius * np.tan(rad_theta / 2) - Parameters - ---------- - ax : matplotlib.axes.Axes, optional - Axis on which to draw the geometry. If None, a new figure and axis are - created. - show : bool, optional - If True, the plot is displayed immediately using plt.show(). - Default is False. - Returns - ------- - matplotlib.axes.Axes - The axis object containing the plot. +@dataclass +class CaseGeometryOptVariables(OptVariablesFrame): + """Optimisiation variables for Trapezoidal Geometry.""" - Notes - ----- - Must be implemented by each specific geometry class. - """ + Ri: OptVariable = ov( + "Ri", + 0, + lower_bound=0, + upper_bound=np.inf, + description="External radius of the TF coil case [m].", + ) + Rk: OptVariable = ov( + "Rk", + 0, + lower_bound=0, + upper_bound=np.inf, + description="Internal radius of the TF coil case [m].", + ) + theta_TF: OptVariable = ov( + "theta_TF", + 0, + lower_bound=0, + upper_bound=360, + description="Toroidal angular span of the TF coil [degrees].", + ) -class TrapezoidalGeometry(CaseGeometry): +class TrapezoidalGeometry(GeometryParameterisation[CaseGeometryOptVariables]): """ Geometry of a Toroidal Field (TF) coil case with trapezoidal cross-section. The coil cross-section has a trapezoidal shape: wider at the outer radius (Ri) and narrower at the inner radius (Rk), reflecting typical TF coil designs - for magnetic and mechanical optimization. + for magnetic and mechanical optimisation. """ + def __init__(self, var_dict: VarDictT | None = None): + variables = CaseGeometryOptVariables() + variables.adjust_variables(var_dict, strict_bounds=False) + super().__init__(variables) + + @property + def rad_theta(self) -> float: + """ + Compute the Toroidal angular span of the TF coil in radians + """ + return np.radians(self.variables.theta_TF.value) + @property def area(self) -> float: """ @@ -255,66 +121,45 @@ def area(self) -> float: Returns ------- - float + : Cross-sectional area [m²]. """ return ( 0.5 - * (self.dx_at_radius(self.Ri) + self.dx_at_radius(self.Rk)) - * (self.Ri - self.Rk) + * (self.variables.Ri.value - self.variables.Rk.value) + * ( + _dx_at_radius(self.variables.Ri.value, self.rad_theta) + + _dx_at_radius(self.variables.Rk.value, self.rad_theta) + ) ) - def build_polygon(self) -> np.ndarray: + def create_shape(self, label: str = "") -> BluemiraWire: """ Construct the (x, r) coordinates of the trapezoidal cross-section polygon. Returns ------- - np.ndarray + : Array of shape (4, 2) representing the corners of the trapezoid. Coordinates are ordered counterclockwise starting from the top-left corner: - [(-dx_outer/2, Ri), (dx_outer/2, Ri), (dx_inner/2, Rk), (-dx_inner/2, Rk)]. - """ - dx_outer = self.dx_at_radius(self.Ri) - dx_inner = self.dx_at_radius(self.Rk) - - return np.array([ - [-dx_outer / 2, self.Ri], - [dx_outer / 2, self.Ri], - [dx_inner / 2, self.Rk], - [-dx_inner / 2, self.Rk], - ]) - - def plot(self, ax=None, *, show=False) -> plt.Axes: - """ - Plot the trapezoidal cross-sectional shape of the TF case. - - Parameters - ---------- - ax : matplotlib.axes.Axes, optional - Axis object on which to draw the geometry. If None, a new figure and axis - are created. - show : bool, optional - If True, the plot is immediately displayed using plt.show(). Default is - False. - - Returns - ------- - matplotlib.axes.Axes - Axis object containing the plotted geometry. - """ - if ax is None: - _, ax = plt.subplots() - poly = self.build_polygon() - poly = np.vstack([poly, poly[0]]) # Close the polygon - ax.plot(poly[:, 0], poly[:, 1], "k-", linewidth=2) - ax.set_aspect("equal") - if show: - plt.show() - return ax + [(-dx_outer, Ri), (dx_outer, Ri), (dx_inner, Rk), (-dx_inner, Rk)]. + """ + dx_outer = _dx_at_radius(self.variables.Ri.value, self.rad_theta) / 2 + dx_inner = _dx_at_radius(self.variables.Rk.value, self.rad_theta) / 2 + + return make_polygon( + [ + [-dx_outer, 0.0, self.variables.Ri.value], + [dx_outer, 0.0, self.variables.Ri.value], + [dx_inner, 0.0, self.variables.Rk.value], + [-dx_inner, 0.0, self.variables.Rk.value], + ], + closed=True, + label=label, + ) -class WedgedGeometry(CaseGeometry): +class WedgedGeometry(GeometryParameterisation[CaseGeometryOptVariables]): """ TF coil case shaped as a sector of an annulus (wedge with arcs). @@ -322,241 +167,145 @@ class WedgedGeometry(CaseGeometry): connected by radial lines, forming a wedge-like shape. """ + def __init__(self, var_dict: VarDictT | None = None): + variables = CaseGeometryOptVariables() + variables.adjust_variables(var_dict, strict_bounds=False) + super().__init__(variables) + + @property + def rad_theta(self) -> float: + """ + Compute the Toroidal angular span of the TF coil in radians + """ + return np.radians(self.variables.theta_TF.value) + def area(self) -> float: """ Compute the cross-sectional area of the wedge geometry. Returns ------- - float + : Cross-sectional area [m²] defined by the wedge between outer radius Ri and inner radius Rk over the toroidal angle theta_TF. """ - return 0.5 * self.rad_theta_TF * (self.Ri**2 - self.Rk**2) + return ( + 0.5 + * self.rad_theta + * (self.variables.Ri.value**2 - self.variables.Rk.value**2) + ) - def build_polygon(self, n_points: int = 50) -> np.ndarray: + def create_shape(self, label: str = "", n_points: int = 50) -> BluemiraWire: """ Build the polygon representing the wedge shape. - The polygon is created by discretizing the outer and inner arcs + The polygon is created by discretising the outer and inner arcs into a series of points connected sequentially. Parameters ---------- - n_points : int, optional + n_points: Number of points to discretize each arc. Default is 50. Returns ------- - np.ndarray + : Array of (x, y) coordinates [m] describing the wedge polygon. """ - theta1 = -self.rad_theta_TF / 2 + theta1 = -self.rad_theta / 2 theta2 = -theta1 angles_outer = np.linspace(theta1, theta2, n_points) angles_inner = np.linspace(theta2, theta1, n_points) arc_outer = np.column_stack(( - self.Ri * np.sin(angles_outer), - self.Ri * np.cos(angles_outer), + self.variables.Ri.value * np.sin(angles_outer), + self.variables.Ri.value * np.cos(angles_outer), )) arc_inner = np.column_stack(( - self.Rk * np.sin(angles_inner), - self.Rk * np.cos(angles_inner), + self.variables.Rk.value * np.sin(angles_inner), + self.variables.Rk.value * np.cos(angles_inner), )) - return np.vstack((arc_outer, arc_inner)) - - def plot(self, ax=None, *, show=False): - """ - Plot the wedge-shaped TF coil case cross-section. - - Parameters - ---------- - ax : matplotlib.axes.Axes, optional - Axis on which to draw the geometry. If None, a new figure and axis are - created. - show : bool, optional - If True, immediately display the plot with plt.show(). Default is False. - - Returns - ------- - matplotlib.axes.Axes - The axis object containing the plot. - """ - if ax is None: - _, ax = plt.subplots() - poly = self.build_polygon() - poly = np.vstack([poly, poly[0]]) # Close the polygon - ax.plot(poly[:, 0], poly[:, 1], "k-", linewidth=2) - ax.set_aspect("equal") - if show: - plt.show() - return ax + return make_polygon(np.vstack((arc_outer, arc_inner)), label=label) -# ------------------------------------------------------------------------------ -# CaseTF Class -# ------------------------------------------------------------------------------ -class BaseCaseTF(CaseGeometry, ABC, metaclass=RegistrableMeta): +class CaseTF(ABC): """ Abstract Base Class for Toroidal Field Coil Case configurations. Defines the universal properties common to all TF case geometries. """ - _registry_ = CASETF_REGISTRY - _name_in_registry_ = None - def __init__( self, Ri: float, # noqa: N803 + theta_TF: float, dy_ps: float, dy_vault: float, - theta_TF: float, mat_case: Material, - WPs: list[WindingPack], # noqa: N803 + wps: list[WindingPack], + geometry: GeometryParameterisation, name: str = "BaseCaseTF", ): """ - Initialize a BaseCaseTF instance. + Initialise a BaseCaseTF instance. Parameters ---------- - Ri : float - External radius at the top of the TF coil case [m]. - dy_ps : float + Ri: + External radius of the TF coil case [m]. + theta_TF: + Toroidal angular span of the TF coil [degrees]. + dy_ps: Radial thickness of the poloidal support region [m]. - dy_vault : float + dy_vault: Radial thickness of the vault support region [m]. - theta_TF : float - Toroidal angular aperture of the coil [degrees]. - mat_case : Material + mat_case: Structural material assigned to the TF coil case. - WPs : list[WindingPack] + wps: List of winding pack objects embedded inside the TF case. - name : str, optional + name: String identifier for the TF coil case instance (default is "BaseCaseTF"). """ - self._name = None - self.name = name - - self._dy_ps = None + self.geometry = geometry + self.geometry.variables.Ri.value = Ri + self.geometry.variables.theta_TF.value = theta_TF self.dy_ps = dy_ps - - self._WPs = None - self.WPs = WPs - - self._mat_case = None - self.mat_case = mat_case - - self._Ri = None - self.Ri = Ri - - self._theta_TF = None - self.theta_TF = theta_TF - - # super().__init__(Ri=Ri, Rk=0, theta_TF=theta_TF) - - self._dy_vault = None self.dy_vault = dy_vault + self.mat_case = mat_case + self.wps = wps + self.name = name + # sets Rk + self.update_dy_vault(self.dy_vault) @property - def name(self) -> str: - """ - Name identifier of the TF case. - - Returns - ------- - str - Human-readable label for the coil case instance. - """ - return self._name - - @name.setter - def name(self, value: str): - """ - Set the name of the TF case. - - Parameters - ---------- - value : str - Case name. - - Raises - ------ - TypeError - If value is not a string. - """ - if not isinstance(value, str): - raise TypeError("name must be a string.") - self._name = value - - @property - def dy_ps(self) -> float: + def rad_theta(self) -> float: """ - Radial thickness of the poloidal support (PS) region [m]. - - Returns - ------- - float - Thickness of the upper structural cap between the TF case wall and the - first winding pack [m]. + Compute the Toroidal angular span of the TF coil in radians """ - return self._dy_ps + return np.radians(self.geometry.variables.theta_TF.value) - @dy_ps.setter - def dy_ps(self, value: float): + def update_dy_vault(self, value: float): """ - Set the thickness of the poloidal support region. + Update the value of the vault support region thickness Parameters ---------- - value : float - Poloidal support thickness [m]. - - Raises - ------ - ValueError - If value is not positive. - """ - if value <= 0: - raise ValueError("dy_ps must be positive.") - self._dy_ps = value - - @property - def dy_vault(self) -> float: - """ - Radial thickness of the vault support region [m]. - - Returns - ------- - float - Thickness of the lower structural region supporting the winding packs [m]. + value: + Vault thickness [m]. """ - return self._dy_vault + self.dy_vault = value + self.geometry.variables.Rk.value = self.R_wp_k[-1] - self.dy_vault - @dy_vault.setter - def dy_vault(self, value: float): + # jm - not used but replaces functionality of original Rk setter + # can't find when (if) it was used originally + def update_Rk(self, value: float): # noqa: N802 """ - Set the thickness of the vault support region. - - Parameters - ---------- - value : float - Vault thickness [m]. - - Raises - ------ - ValueError - If value is not positive. + Set or update the internal (innermost) radius of the TF case. """ - if value <= 0: - raise ValueError("dy_vault must be positive.") - self._dy_vault = value - - self.Rk = self.R_wp_k[-1] - self._dy_vault + self.geometry.variables.Rk.value = value + self.dy_vault = self.R_wp_k[-1] - self.geometry.variables.Rk.value @property @abstractmethod @@ -566,7 +315,7 @@ def dx_vault(self): Returns ------- - float + : Average length of the vault in the toroidal direction [m]. """ @@ -577,7 +326,7 @@ def mat_case(self) -> Material: Returns ------- - Material + : Material object providing mechanical and thermal properties. """ return self._mat_case @@ -589,7 +338,7 @@ def mat_case(self, value: Material): Parameters ---------- - value : Material + value: Material object. Raises @@ -601,25 +350,25 @@ def mat_case(self, value: Material): self._mat_case = value @property - def WPs(self) -> list[WindingPack]: # noqa: N802 + def wps(self) -> list[WindingPack]: """ List of winding pack (WP) objects embedded inside the TF case. Returns ------- - list of WindingPack + : Winding pack instances composing the internal coil layout. """ - return self._WPs + return self._wps - @WPs.setter - def WPs(self, value: list[WindingPack]): # noqa: N802 + @wps.setter + def wps(self, value: list[WindingPack]): """ Set the winding pack objects list. Parameters ---------- - value : list[WindingPack] + value: List containing only WindingPack objects. Raises @@ -628,29 +377,31 @@ def WPs(self, value: list[WindingPack]): # noqa: N802 If value is not a list of WindingPack instances. """ if not isinstance(value, list): - raise TypeError("WPs must be a list of WindingPack objects.") + raise TypeError("wps must be a list of WindingPack objects.") if not all(isinstance(wp, WindingPack) for wp in value): - raise TypeError("All elements of WPs must be WindingPack instances.") - self._WPs = value + raise TypeError("All elements of wps must be WindingPack instances.") + self._wps = value # fix dy_vault (this will recalculate Rk) - if hasattr(self, "dy_vault"): - self.dy_vault = self.dy_vault + self.update_dy_vault(self.dy_vault) @property def dx_i(self): """Toroidal length of the coil case at its maximum radial position [m]""" - return 2 * self.Ri * np.tan(self._rad_theta_TF / 2) + return _dx_at_radius(self.geometry.variables.Ri.value, self.rad_theta) @property def dx_ps(self): """Average toroidal length of the ps plate [m]""" - return (self.Ri + (self.Ri - self.dy_ps)) * np.tan(self._rad_theta_TF / 2) + return ( + self.geometry.variables.Ri.value + + (self.geometry.variables.Ri.value - self.dy_ps) + ) * np.tan(self.rad_theta / 2) @property - def n_conductors(self): + def n_conductors(self) -> int: """Total number of conductors in the winding pack.""" - return sum(w.n_conductors for w in self.WPs) + return sum(w.n_conductors for w in self.wps) @property def dy_wp_i(self) -> np.ndarray: @@ -659,11 +410,11 @@ def dy_wp_i(self) -> np.ndarray: Returns ------- - np.ndarray + : Array containing the radial thickness [m] of each Winding Pack. - Each element corresponds to one WP in the self.WPs list. + Each element corresponds to one WP in the self.wps list. """ - return np.array([wp.dy for wp in self.WPs]) + return np.array([wp.dy for wp in self.wps]) @property def dy_wp_tot(self) -> float: @@ -672,7 +423,7 @@ def dy_wp_tot(self) -> float: Returns ------- - float + : Total radial thickness [m] summed over all winding packs. """ return sum(self.dy_wp_i) @@ -684,18 +435,18 @@ def R_wp_i(self) -> np.ndarray: # noqa: N802 Returns ------- - np.ndarray + : Array of radial positions [m] corresponding to the outer edge of each WP. """ dy_wp_cumsum = np.cumsum(np.concatenate(([0.0], self.dy_wp_i))) - result_initial = self.Ri - self.dy_ps + result_initial = self.geometry.variables.Ri.value - self.dy_ps if len(dy_wp_cumsum) == 1: result = np.array([result_initial]) else: result = result_initial - dy_wp_cumsum[:-1] - if len(result) != len(self.WPs): - bluemira_error(f"Mismatch: {len(result)} R_wp_i vs {len(self.WPs)} WPs!") + if len(result) != len(self.wps): + bluemira_error(f"Mismatch: {len(result)} R_wp_i vs {len(self.wps)} wps!") return result @@ -706,87 +457,53 @@ def R_wp_k(self): # noqa: N802 Returns ------- - np.ndarray + : Array of radial positions [m] corresponding to the outer edge of each winding pack. """ return self.R_wp_i - self.dy_wp_i - @property - def Rk(self) -> float: # noqa: N802 - """ - Internal (innermost) radius of the TF case at the top [m]. - - Returns - ------- - float - Inner radius measured from the machine center to the case outer wall [m]. - """ - return self._Rk - - @Rk.setter - def Rk(self, value: float): # noqa: N802 - """ - Set the internal (innermost) radius of the TF case. - - Parameters - ---------- - value : float - Outer radius [m]. Must be a strictly positive number. - - Raises - ------ - ValueError - If the provided radius is not positive. - """ - if value < 0: - raise ValueError("Rk must be positive.") - self._Rk = value - - self._dy_vault = self.R_wp_k[-1] - self._Rk - - def plot(self, ax=None, *, show: bool = False, homogenized: bool = False): + def plot( + self, + ax: plt.Axes | None = None, + *, + show: bool = False, + homogenised: bool = False, + ) -> plt.Axes: """ Schematic plot of the TF case cross-section including winding packs. Parameters ---------- - ax : matplotlib.axes.Axes, optional + ax: Axis on which to draw the figure. If `None`, a new figure and axis will be created. - show : bool, optional + show: If `True`, displays the plot immediately using `plt.show()`. Default is `False`. - homogenized : bool, optional - If `True`, plots winding packs as homogenized blocks. - If `False`, plots individual conductors inside WPs. + homogenised: + If `True`, plots winding packs as homogenised blocks. + If `False`, plots individual conductors inside wps. Default is `False`. Returns ------- - matplotlib.axes.Axes + : The axis object containing the rendered plot. """ if ax is None: _, ax = plt.subplots() ax.set_aspect("equal", adjustable="box") - # -------------------------------------- - # Plot external case boundary (delegate) - # -------------------------------------- - super().plot(ax=ax, show=False) + self.geometry.plot(ax=ax) - # -------------------------------------- # Plot winding packs - # -------------------------------------- - for i, wp in enumerate(self.WPs): + for i, wp in enumerate(self.wps): xc_wp = 0.0 yc_wp = self.R_wp_i[i] - wp.dy / 2 - ax = wp.plot(xc=xc_wp, yc=yc_wp, ax=ax, homogenized=homogenized) + ax = wp.plot(xc=xc_wp, yc=yc_wp, ax=ax, homogenised=homogenised) - # -------------------------------------- - # Finalize plot - # -------------------------------------- + # Finalise plot ax.set_xlabel("Toroidal direction [m]") ax.set_ylabel("Radial direction [m]") ax.set_title(f"TF Case Cross Section: {self.name}") @@ -797,40 +514,39 @@ def plot(self, ax=None, *, show: bool = False, homogenized: bool = False): return ax @property - def area_case_jacket(self): + def area_case_jacket(self) -> float: """ Area of the case jacket (excluding winding pack regions). Returns ------- - float Case jacket area [m²], computed as total area minus total WP area. """ - return self.area - self.area_wps + return self.geometry.area - self.area_wps @property - def area_wps(self): + def area_wps(self) -> float: """ Total area occupied by all winding packs. Returns ------- - float + : Combined area of the winding packs [m²]. """ - return np.sum([w.area for w in self.WPs]) + return np.sum([w.area for w in self.wps]) @property - def area_wps_jacket(self): + def area_wps_jacket(self) -> float: """ Total jacket area of all winding packs. Returns ------- - float - Combined area of conductor jackets in all WPs [m²]. + : + Combined area of conductor jackets in all wps [m²]. """ - return np.sum([w.jacket_area for w in self.WPs]) + return np.sum([w.jacket_area for w in self.wps]) @property def area_jacket_total(self) -> float: @@ -839,11 +555,11 @@ def area_jacket_total(self) -> float: - The case jacket area (structural material surrounding the winding packs). - The conductor jackets area (jackets enclosing the individual conductors - inside the WPs). + inside the wps). Returns ------- - float + : Combined area of the case structure and the conductor jackets [m²]. Notes @@ -867,15 +583,15 @@ def rearrange_conductors_in_wp( Parameters ---------- - n_conductors : int + n_conductors: Total number of conductors to distribute. - wp_reduction_factor : float - Fractional reduction of available toroidal space for WPs. - min_gap_x : float + wp_reduction_factor: + Fractional reduction of available toroidal space for wps. + min_gap_x: Minimum gap between the WP and the case boundary in toroidal direction [m]. - n_layers_reduction : int + n_layers_reduction: Number of layers to remove after each WP. - layout : str, optional + layout: Layout strategy ("auto", "layer", "pancake"). """ @@ -895,15 +611,15 @@ def enforce_wp_layout_rules( Parameters ---------- - n_conductors : int + n_conductors: Number of conductors to allocate. - dx_WP : float + dx_WP: Available toroidal width for the winding pack [m]. - dx_cond : float + dx_cond: Toroidal width of a single conductor [m]. - dy_cond : float + dy_cond: Radial height of a single conductor [m]. - layout : str + layout: Layout type: - "auto" : no constraints - "layer" : enforce even number of turns (ny % 2 == 0) @@ -911,8 +627,10 @@ def enforce_wp_layout_rules( Returns ------- - tuple[int, int] - n_layers_max (nx), n_turns_max (ny) + : + n_layers_max (nx) + : + n_turns_max (ny) Raises ------ @@ -947,165 +665,88 @@ def enforce_wp_layout_rules( return n_layers_max, n_turns_max - @abstractmethod - def optimize_vault_radial_thickness( - self, - pm: float, - fz: float, - T: float, # noqa: N803 - B: float, - allowable_sigma: float, - bounds: np.ndarray = None, - ): - """ - Abstract method to optimize the radial thickness of the vault support region. - - Parameters - ---------- - pm : float - Radial magnetic pressure [Pa]. - fz : float - Axial electromagnetic force [N]. - T : float - Operating temperature [K]. - B : float - Magnetic field strength [T]. - allowable_sigma : float - Allowable maximum stress [Pa]. - bounds : np.ndarray, optional - Optimization bounds for vault thickness [m]. - """ - - def to_dict(self) -> dict: + def to_dict(self) -> dict[str, float | str | list[dict[str, float | str | Any]]]: """ - Serialize the BaseCaseTF instance into a dictionary. + Serialise the BaseCaseTF instance into a dictionary. Returns ------- dict - Serialized data representing the TF case, including geometry and material + Serialised data representing the TF case, including geometry and material information. """ return { - "name_in_registry": getattr( - self, "_name_in_registry_", self.__class__.__name__ - ), "name": self.name, - "Ri": self.Ri, + "Ri": self.geometry.variables.Ri.value, "dy_ps": self.dy_ps, "dy_vault": self.dy_vault, - "theta_TF": self.theta_TF, - "mat_case": self.mat_case.name, # Assume Material has 'name' attribute - "WPs": [wp.to_dict() for wp in self.WPs], - # Assume each WindingPack implements to_dict() + "theta_TF": self.geometry.variables.theta_TF.value, + "mat_case": self.mat_case, + "wps": [wp.to_dict() for wp in self.wps], } - @classmethod - def from_dict(cls, case_dict: dict, name: str | None = None) -> "BaseCaseTF": - """ - Deserialize a BaseCaseTF instance from a dictionary. - - Parameters - ---------- - case_dict : dict - Dictionary containing serialized TF case data. - name : str, optional - Optional name override for the new instance. - - Returns - ------- - BaseCaseTF - Reconstructed TF case instance. - - Raises - ------ - ValueError - If the 'name_in_registry' field does not match this class. - """ - name_in_registry = case_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - WPs = [create_wp_from_dict(wp_dict) for wp_dict in case_dict["WPs"]] # noqa:N806 - - return cls( - Ri=case_dict["Ri"], - dy_ps=case_dict["dy_ps"], - dy_vault=case_dict["dy_vault"], - theta_TF=case_dict["theta_TF"], - mat_case=case_dict["mat_case"], - WPs=WPs, - name=name or case_dict.get("name"), - ) - def __str__(self) -> str: """ Generate a human-readable summary of the TF case. Returns ------- - str - Multiline string summarizing key properties of the TF case. + : + Multiline string summarising key properties of the TF case. """ return ( f"CaseTF '{self.name}'\n" - f" - Ri: {self.Ri:.3f} m\n" - f" - Rk: {self.Rk:.3f} m\n" + f" - Ri: {self.geometry.variables.Ri.value:.3f} m\n" + f" - Rk: {self.geometry.variables.Rk.value:.3f} m\n" f" - dy_ps: {self.dy_ps:.3f} m\n" f" - dy_vault: {self.dy_vault:.3f} m\n" - f" - theta_TF: {self.theta_TF:.2f}°\n" + f" - theta_TF: {self.geometry.variables.theta_TF.value:.2f}°\n" f" - Material: {self.mat_case.name}\n" - f" - Winding Packs: {len(self.WPs)} packs\n" + f" - Winding Packs: {len(self.wps)} packs\n" ) -class TrapezoidalCaseTF(BaseCaseTF, TrapezoidalGeometry): +class TrapezoidalCaseTF(CaseTF): """ Toroidal Field Coil Case with Trapezoidal Geometry. Note: this class considers a set of Winding Pack with the same conductor (instance). """ - _registry_ = CASETF_REGISTRY - _name_in_registry_ = "TrapezoidalCaseTF" - def __init__( self, Ri: float, # noqa: N803 + theta_TF: float, dy_ps: float, dy_vault: float, - theta_TF: float, mat_case: Material, - WPs: list[WindingPack], # noqa: N803 + wps: list[WindingPack], name: str = "TrapezoidalCaseTF", ): - self._check_WPs(WPs) + self._check_wps(wps) + geom = TrapezoidalGeometry() super().__init__( Ri=Ri, + theta_TF=theta_TF, dy_ps=dy_ps, dy_vault=dy_vault, - theta_TF=theta_TF, mat_case=mat_case, - WPs=WPs, + wps=wps, name=name, + geometry=geom, ) - def _check_WPs( # noqa: PLR6301, N802 + def _check_wps( # noqa: PLR6301 self, - WPs: list[WindingPack], # noqa:N803 + wps: list[WindingPack], ): """ - Validate that the provided winding packs (WPs) are non-empty and share the + Validate that the provided winding packs (wps) are non-empty and share the same conductor. Parameters ---------- - WPs : list of WindingPack + wps: List of winding pack objects to validate. Raises @@ -1115,15 +756,15 @@ def _check_WPs( # noqa: PLR6301, N802 ValueError If winding packs have different conductor instances. """ - if not WPs: + if not wps: raise ValueError("At least one non-empty winding pack must be provided.") - first_conductor = WPs[0].conductor - for i, wp in enumerate(WPs[1:], start=1): + first_conductor = wps[0].conductor + for i, wp in enumerate(wps[1:], start=1): if wp.conductor is not first_conductor: bluemira_warn( f"[Winding pack at index {i} uses a different conductor object " - f"than the first one. This module requires all WPs to " + f"than the first one. This module requires all wps to " f"share the same conductor instance." f"Please verify the inputs or unify the conductor assignment." ) @@ -1135,10 +776,12 @@ def dx_vault(self): Returns ------- - float + : Average length of the vault in the toroidal direction [m]. """ - return (self.R_wp_k[-1] + self.Rk) * np.tan(self.rad_theta_TF / 2) + return (self.R_wp_k[-1] + self.geometry.variables.Rk.value) * np.tan( + self.rad_theta / 2 + ) def Kx_ps(self, op_cond: OperationalConditions): # noqa: N802 """ @@ -1176,31 +819,30 @@ def Kx_lat(self, op_cond: OperationalConditions): # noqa: N802 Array of radial stiffness values for each lateral segment [Pa]. """ dx_lat = np.array([ - (self.R_wp_i[i] + self.R_wp_k[i]) / 2 * np.tan(self.rad_theta_TF / 2) - - w.dx / 2 - for i, w in enumerate(self.WPs) + (self.R_wp_i[i] + self.R_wp_k[i]) / 2 * np.tan(self.rad_theta / 2) - w.dx / 2 + for i, w in enumerate(self.wps) ]) - dy_lat = np.array([w.dy for w in self.WPs]) + dy_lat = np.array([w.dy for w in self.wps]) return self.mat_case.youngs_modulus(op_cond) * dy_lat / dx_lat - def Kx_vault(self, op_cond: OperationalConditions): # noqa: N802 + def Kx_vault(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the equivalent radial stiffness of the vault region. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Equivalent radial stiffness of the vault [Pa]. """ return self.mat_case.youngs_modulus(op_cond) * self.dy_vault / self.dx_vault - def Kx(self, op_cond: OperationalConditions): # noqa: N802 + def Kx(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the total equivalent radial stiffness of the entire case structure. @@ -1211,24 +853,26 @@ def Kx(self, op_cond: OperationalConditions): # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Total equivalent radial stiffness of the TF case [Pa]. """ + # TODO possible reason for floating point difference + kx_lat = self.Kx_lat(op_cond) temp = [ - serie_k([ - self.Kx_lat(op_cond)[i], + reciprocal_summation([ + kx_lat[i], w.Kx(op_cond), - self.Kx_lat(op_cond)[i], + kx_lat[i], ]) - for i, w in enumerate(self.WPs) + for i, w in enumerate(self.wps) ] - return parall_k([self.Kx_ps(op_cond), self.Kx_vault(op_cond), *temp]) + return summation([self.Kx_ps(op_cond), self.Kx_vault(op_cond), *temp]) def Ky_ps(self, op_cond: OperationalConditions): # noqa: N802 """ @@ -1242,7 +886,7 @@ def Ky_ps(self, op_cond: OperationalConditions): # noqa: N802 Returns ------- - float + : Equivalent toroidal stiffness of the PS region [Pa]. """ return self.mat_case.youngs_modulus(op_cond) * self.dx_ps / self.dy_ps @@ -1266,11 +910,11 @@ def Ky_lat(self, op_cond: OperationalConditions): # noqa: N802 Array of toroidal stiffness values for each lateral segment [Pa]. """ dx_lat = np.array([ - (self.R_wp_i[i] + self.R_wp_k[i]) / 2 * np.tan(self._rad_theta_TF / 2) + (self.R_wp_i[i] + self.R_wp_k[i]) / 2 * np.tan(self._rad_theta / 2) - w.dx / 2 - for i, w in enumerate(self.WPs) + for i, w in enumerate(self.wps) ]) - dy_lat = np.array([w.dy for w in self.WPs]) + dy_lat = np.array([w.dy for w in self.wps]) return self.mat_case.youngs_modulus(op_cond) * dx_lat / dy_lat def Ky_vault(self, op_cond: OperationalConditions): # noqa: N802 @@ -1290,7 +934,7 @@ def Ky_vault(self, op_cond: OperationalConditions): # noqa: N802 """ return self.mat_case.youngs_modulus(op_cond) * self.dx_vault / self.dy_vault - def Ky(self, op_cond: OperationalConditions): # noqa: N802 + def Ky(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Compute the total equivalent toroidal stiffness of the entire case structure. @@ -1301,24 +945,26 @@ def Ky(self, op_cond: OperationalConditions): # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Total equivalent toroidal stiffness of the TF case [Pa]. """ + # TODO possible reason for floating point difference + ky_lat = self.Ky_lat(op_cond) temp = [ - parall_k([ - self.Ky_lat(op_cond)[i], + summation([ + ky_lat[i], w.Ky(op_cond), - self.Ky_lat(op_cond)[i], + ky_lat[i], ]) - for i, w in enumerate(self.WPs) + for i, w in enumerate(self.wps) ] - return serie_k([self.Ky_ps(op_cond), self.Ky_vault(op_cond), *temp]) + return reciprocal_summation([self.Ky_ps(op_cond), self.Ky_vault(op_cond), *temp]) def rearrange_conductors_in_wp( self, @@ -1329,20 +975,20 @@ def rearrange_conductors_in_wp( layout: str = "auto", ): """ - Rearrange the total number of conductors into winding packs (WPs) + Rearrange the total number of conductors into winding packs within the TF coil case geometry using enforce_wp_layout_rules. Parameters ---------- - n_conductors : int + n_conductors: Total number of conductors to be allocated. - wp_reduction_factor : float - Fractional reduction of the total available toroidal space for WPs. - min_gap_x : float + wp_reduction_factor: + Fractional reduction of the total available toroidal space for winding packs. + min_gap_x: Minimum allowable toroidal gap between WP and boundary [m]. - n_layers_reduction : int + n_layers_reduction: Number of horizontal layers to reduce after each WP. - layout : str, optional + layout: Layout type ("auto", "layer", "pancake"). Raises @@ -1351,7 +997,7 @@ def rearrange_conductors_in_wp( If there is not enough space to allocate all the conductors. """ debug_msg = ["Method rearrange_conductors_in_wp"] - conductor = self.WPs[0].conductor + conductor = self.wps[0].conductor R_wp_i = self.R_wp_i[0] # noqa: N806 dx_WP = self.dx_i * wp_reduction_factor # noqa: N806 @@ -1365,22 +1011,18 @@ def rearrange_conductors_in_wp( f"n_conductors = {n_conductors}", ]) - WPs = [] # noqa: N806 + wps = [] # number of conductors to be allocated remaining_conductors = n_conductors - # maximum number of winding packs in WPs + # maximum number of winding packs in wps i_max = 50 - i = 0 - while i < i_max and remaining_conductors > 0: - i += 1 - + n_layers_max = 0 + for i in range(i_max): # maximum number of turns on the considered WP - if i == 1: + if i == 0: n_layers_max = math.floor(dx_WP / conductor.dx) if layout == "pancake": - n_layers_max = math.floor(dx_WP / conductor.dx / 2.0) * 2 - if n_layers_max == 0: - n_layers_max = 2 + n_layers_max = (math.floor(dx_WP / conductor.dx / 2) * 2) or 2 else: n_layers_max -= n_layers_reduction @@ -1391,17 +1033,17 @@ def rearrange_conductors_in_wp( ) if n_layers_max >= remaining_conductors: - WPs.append( + wps.append( WindingPack(conductor=conductor, nx=remaining_conductors, ny=1) ) remaining_conductors = 0 else: dx_WP = n_layers_max * conductor.dx # noqa: N806 - gap_0 = R_wp_i * np.tan(self.rad_theta_TF / 2) - dx_WP / 2 + gap_0 = R_wp_i * np.tan(self.rad_theta / 2) - dx_WP / 2 gap_1 = min_gap_x - max_dy = (gap_0 - gap_1) / np.tan(self.rad_theta_TF / 2) + max_dy = (gap_0 - gap_1) / np.tan(self.rad_theta / 2) n_turns_max = min( int(np.floor(max_dy / conductor.dy)), int(np.ceil(remaining_conductors / n_layers_max)), @@ -1422,34 +1064,35 @@ def rearrange_conductors_in_wp( if n_layers_max * n_turns_max > remaining_conductors: n_turns_max -= 1 - WPs.append( + wps.append( WindingPack(conductor=conductor, nx=n_layers_max, ny=n_turns_max) ) remaining_conductors -= n_layers_max * n_turns_max - WPs.append( + wps.append( WindingPack(conductor=conductor, nx=remaining_conductors, ny=1) ) remaining_conductors = 0 else: - WPs.append( + wps.append( WindingPack(conductor=conductor, nx=n_layers_max, ny=n_turns_max) ) remaining_conductors -= n_layers_max * n_turns_max - if remaining_conductors < 0: - bluemira_warn( - f"{abs(remaining_conductors)}/{n_layers_max * n_turns_max}" - f"have been added to complete the last winding pack (nx" - f"={n_layers_max}, ny={n_turns_max})." - ) - R_wp_i -= n_turns_max * conductor.dy # noqa: N806 debug_msg.append( f"n_layers_max: {n_layers_max}, n_turns_max: {n_turns_max}" ) + if remaining_conductors <= 0: + if remaining_conductors < 0: + bluemira_warn( + f"{abs(remaining_conductors)}/{n_layers_max * n_turns_max}" + f"have been added to complete the last winding pack (nx" + f"={n_layers_max}, ny={n_turns_max})." + ) + break bluemira_debug("\n".join(debug_msg)) - self.WPs = WPs + self.wps = wps # just a final check if self.n_conductors != n_conductors: @@ -1461,7 +1104,9 @@ def rearrange_conductors_in_wp( bluemira_error(msg) raise ValueError(msg) - def _tresca_stress(self, pm: float, fz: float, op_cond: OperationalConditions): + def _tresca_stress( + self, pm: float, fz: float, op_cond: OperationalConditions + ) -> float: """ Estimate the maximum principal (Tresca) stress on the inner case of the TF coil. @@ -1475,24 +1120,26 @@ def _tresca_stress(self, pm: float, fz: float, op_cond: OperationalConditions): Parameters ---------- - pm : float + pm: Radial magnetic pressure acting on the case [Pa]. - fz : float + fz: Vertical force acting on the inner leg of the case [N]. - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Estimated maximum stress [Pa] acting on the case nose (hoop + vertical contribution). """ # The maximum principal stress acting on the case nose is the compressive # hoop stress generated in the equivalent shell from the magnetic pressure. From # the Shell theory, for an isotropic continuous shell with a thickness ratio: - beta = self.Rk / (self.Rk + self.dy_vault) + beta = self.geometry.variables.Rk.value / ( + self.geometry.variables.Rk.value + self.dy_vault + ) # the maximum hoop stress, corrected to account for the presence of the WP, is # placed at the innermost radius of the case as: sigma_theta = ( @@ -1508,62 +1155,6 @@ def _tresca_stress(self, pm: float, fz: float, op_cond: OperationalConditions): sigma_z = fz / (self.area_case_jacket + self.area_wps_jacket) return sigma_theta + sigma_z - def optimize_vault_radial_thickness( - self, - pm: float, - fz: float, - op_cond, - allowable_sigma: float, - bounds: np.array = None, - ): - """ - Optimize the vault radial thickness of the case - - Parameters - ---------- - pm : - The magnetic pressure applied along the radial direction (Pa). - f_z : - The force applied in the z direction, perpendicular to the case - cross-section (N). - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material properties. - allowable_sigma : - The allowable stress (Pa) for the jacket material. - bounds : - Optional bounds for the jacket thickness optimization (default is None). - - Returns - ------- - The result of the optimization process containing information about the - optimal vault thickness. - - Raises - ------ - ValueError - If the optimization process did not converge. - """ - method = None - if bounds is not None: - method = "bounded" - - result = minimize_scalar( - fun=self._sigma_difference, - args=(pm, fz, op_cond, allowable_sigma), - bounds=bounds, - method=method, - options={"xatol": 1e-4}, - ) - - if not result.success: - raise ValueError("dy_vault optimization did not converge.") - self.dy_vault = result.x - # print(f"Optimal dy_vault: {self.dy_vault}") - # print(f"Tresca sigma: {self._tresca_stress(pm, fz, T=T, B=B) / 1e6} MPa") - - return result - def _sigma_difference( self, dy_vault: float, @@ -1571,314 +1162,40 @@ def _sigma_difference( fz: float, op_cond: OperationalConditions, allowable_sigma: float, - ): + ) -> float: """ - Fitness function for the optimization problem. It calculates the absolute + Fitness function for the optimisation problem. It calculates the absolute difference between the Tresca stress and the allowable stress. Parameters ---------- - dy_vault : + dy_vault: The thickness of the vault in the direction perpendicular to the applied pressure(m). - pm : + pm: The magnetic pressure applied along the radial direction (Pa). - fz : + fz: The force applied in the z direction, perpendicular to the case cross-section (N). - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material properties. - allowable_sigma : + allowable_sigma: The allowable stress (Pa) for the vault material. Returns ------- + : The absolute difference between the calculated Tresca stress and the allowable stress (Pa). Notes ----- - This function modifies the case's vault thickness - using the value provided in jacket_thickness. + This function modifies the case's vault thickness + using the value provided in jacket_thickness. """ self.dy_vault = dy_vault sigma = self._tresca_stress(pm, fz, op_cond) # bluemira_print(f"sigma: {sigma}, allowable_sigma: {allowable_sigma}, # diff: {sigma - allowable_sigma}") return abs(sigma - allowable_sigma) - - def optimize_jacket_and_vault( - self, - pm: float, - fz: float, - op_cond: OperationalConditions, - allowable_sigma: float, - bounds_cond_jacket: np.array = None, - bounds_dy_vault: np.array = None, - layout: str = "auto", - wp_reduction_factor: float = 0.8, - min_gap_x: float = 0.05, - n_layers_reduction: int = 4, - max_niter: int = 10, - eps: float = 1e-8, - n_conds: int | None = None, - ): - """ - Jointly optimize the conductor jacket and case vault thickness - under electromagnetic loading constraints. - - This method performs an iterative optimization of: - - The cross-sectional area of the conductor jacket. - - The vault radial thickness of the TF coil casing. - - The optimization loop continues until the relative change in - jacket area and vault thickness drops below the specified - convergence threshold `eps`, or `max_niter` is reached. - - Parameters - ---------- - pm : float - Radial magnetic pressure on the conductor [Pa]. - fz : float - Axial electromagnetic force on the winding pack [N]. - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material properties. - allowable_sigma : float - Maximum allowable stress for structural material [Pa]. - bounds_cond_jacket : np.ndarray, optional - Min/max bounds for conductor jacket area optimization [m²]. - bounds_dy_vault : np.ndarray, optional - Min/max bounds for the case vault thickness optimization [m]. - layout : str, optional - Cable layout strategy; "auto" or predefined layout name. - wp_reduction_factor : float, optional - Reduction factor applied to WP footprint during conductor rearrangement. - min_gap_x : float, optional - Minimum spacing between adjacent conductors [m]. - n_layers_reduction : int, optional - Number of conductor layers to remove when reducing WP height. - max_niter : int, optional - Maximum number of optimization iterations. - eps : float, optional - Convergence threshold for the combined optimization loop. - n_conds : int, optional - Target total number of conductors in the winding pack. If None, the self - number of conductors is used. - - Notes - ----- - The function modifies the internal state of `conductor` and `self.dy_vault`. - """ - debug_msg = ["Method optimize_jacket_and_vault"] - - # Initialize convergence array - self._convergence_array = [] - - if n_conds is None: - n_conds = self.n_conductors - - conductor = self.WPs[0].conductor - - self._check_WPs(self.WPs) - - i = 0 - err_conductor_area_jacket = 10000 * eps - err_dy_vault = 10000 * eps - tot_err = err_dy_vault + err_conductor_area_jacket - - self._convergence_array.append([ - i, - conductor.dy_jacket, - self.dy_vault, - err_conductor_area_jacket, - err_dy_vault, - self.dy_wp_tot, - self.Ri - self.Rk, - ]) - - damping_factor = 0.3 - - while i < max_niter and tot_err > eps: - i += 1 - debug_msg.append(f"Internal optimazion - iteration {i}") - - # Store current values - cond_dx_jacket0 = conductor.dx_jacket - case_dy_vault0 = self.dy_vault - - debug_msg.append( - f"before optimization: conductor jacket area = {conductor.area_jacket}" - ) - cond_area_jacket0 = conductor.area_jacket - t_z_cable_jacket = ( - fz - * self.area_wps_jacket - / (self.area_case_jacket + self.area_wps_jacket) - / self.n_conductors - ) - conductor.optimize_jacket_conductor( - pm, t_z_cable_jacket, op_cond, allowable_sigma, bounds_cond_jacket - ) - debug_msg.extend([ - f"t_z_cable_jacket: {t_z_cable_jacket}", - f"after optimization: conductor jacket area = {conductor.area_jacket}", - ]) - - conductor.dx_jacket = ( - 1 - damping_factor - ) * cond_dx_jacket0 + damping_factor * conductor.dx_jacket - - err_conductor_area_jacket = ( - abs(conductor.area_jacket - cond_area_jacket0) / cond_area_jacket0 - ) - - self.rearrange_conductors_in_wp( - n_conds, - wp_reduction_factor, - min_gap_x, - n_layers_reduction, - layout=layout, - ) - - debug_msg.append(f"before optimization: case dy_vault = {self.dy_vault}") - self.optimize_vault_radial_thickness( - pm=pm, - fz=fz, - op_cond=op_cond, - allowable_sigma=allowable_sigma, - bounds=bounds_dy_vault, - ) - - self.dy_vault = ( - 1 - damping_factor - ) * case_dy_vault0 + damping_factor * self.dy_vault - - delta_case_dy_vault = abs(self.dy_vault - case_dy_vault0) - err_dy_vault = delta_case_dy_vault / self.dy_vault - tot_err = err_dy_vault + err_conductor_area_jacket - - debug_msg.append( - f"after optimization: case dy_vault = {self.dy_vault}\n" - f"err_dy_jacket = {err_conductor_area_jacket}\n " - f"err_dy_vault = {err_dy_vault}\n " - f"tot_err = {tot_err}" - ) - - # Store iteration results in convergence array - self._convergence_array.append([ - i, - conductor.dy_jacket, - self.dy_vault, - err_conductor_area_jacket, - err_dy_vault, - self.dy_wp_tot, - self.Ri - self.Rk, - ]) - - # final check - if i < max_niter: - bluemira_print( - f"Optimization of jacket and vault reached after " - f"{i} iterations. Total error: {tot_err} < {eps}." - ) - - ax = self.plot(show=False, homogenized=False) - ax.set_title("Case design after optimization") - plt.show() - - else: - bluemira_warn( - f"Maximum number of optimization iterations {max_niter} " - f"reached. A total of {tot_err} > {eps} has been obtained." - ) - - def plot_convergence(self): - """ - Plot the evolution of thicknesses and error values over optimization iterations. - - Raises - ------ - RuntimeError - If no convergence data available - """ - if not hasattr(self, "_convergence_array") or not self._convergence_array: - raise RuntimeError("No convergence data available. Run optimization first.") - - convergence_data = np.array(self._convergence_array) - - iterations = convergence_data[:, 0] - dy_jacket = convergence_data[:, 1] - dy_vault = convergence_data[:, 2] - err_dy_jacket = convergence_data[:, 3] - err_dy_vault = convergence_data[:, 4] - dy_wp_tot = convergence_data[:, 5] - Ri_minus_Rk = convergence_data[:, 6] # noqa: N806 - - _, axs = plt.subplots(2, 1, figsize=(10, 10), sharex=True) - - # Top subplot: Thicknesses - axs[0].plot(iterations, dy_jacket, marker="o", label="dy_jacket [m]") - axs[0].plot(iterations, dy_vault, marker="s", label="dy_vault [m]") - axs[0].plot(iterations, dy_wp_tot, marker="^", label="dy_wp_tot [m]") - axs[0].plot(iterations, Ri_minus_Rk, marker="v", label="Ri - Rk [m]") - axs[0].set_ylabel("Thickness [m]") - axs[0].set_title("Evolution of Jacket, Vault, and WP Thicknesses") - axs[0].legend() - axs[0].grid(visible=True) - - # Bottom subplot: Errors - axs[1].plot(iterations, err_dy_jacket, marker="o", label="err_dy_jacket") - axs[1].plot(iterations, err_dy_vault, marker="s", label="err_dy_vault") - axs[1].set_ylabel("Relative Error") - axs[1].set_xlabel("Iteration") - axs[1].set_title("Evolution of Errors during Optimization") - axs[1].set_yscale("log") # Log scale for better visibility if needed - axs[1].legend() - axs[1].grid(visible=True) - - plt.tight_layout() - plt.show() - - -def create_case_tf_from_dict( - case_dict: dict, - name: str | None = None, -) -> BaseCaseTF: - """ - Factory function to create a CaseTF (or subclass) from a serialized dictionary. - - Parameters - ---------- - case_dict : dict - Serialized case dictionary, must include 'name_in_registry' field. - name : str, optional - Name to assign to the created case. If None, uses the name in the dictionary. - - Returns - ------- - BaseCaseTF - A fully instantiated CaseTF (or subclass) object. - - Raises - ------ - ValueError - If no class is registered with the given name_in_registry. - """ - name_in_registry = case_dict.get("name_in_registry") - if name_in_registry is None: - raise ValueError("CaseTF dictionary must include 'name_in_registry' field.") - - case_cls = CASETF_REGISTRY.get(name_in_registry) - if case_cls is None: - available = list(CASETF_REGISTRY.keys()) - raise ValueError( - f"No registered CaseTF class with name_in_registry '{name_in_registry}'. " - f"Available: {available}" - ) - - return case_cls.from_dict( - name=name, - case_dict=case_dict, - ) diff --git a/bluemira/magnets/conductor.py b/bluemira/magnets/conductor.py index 6a7cfa5f70..02277f93af 100644 --- a/bluemira/magnets/conductor.py +++ b/bluemira/magnets/conductor.py @@ -6,45 +6,27 @@ """Conductor class""" -from typing import Any +from __future__ import annotations + +from typing import TYPE_CHECKING, Any import matplotlib.pyplot as plt import numpy as np -from matproplib import OperationalConditions -from matproplib.material import Material -from scipy.optimize import minimize_scalar - -from bluemira.base.look_and_feel import bluemira_debug -from bluemira.magnets.cable import ABCCable, create_cable_from_dict -from bluemira.magnets.registry import RegistrableMeta -from bluemira.magnets.utils import ( - parall_k, - parall_r, - serie_k, - serie_r, -) - -# ------------------------------------------------------------------------------ -# Global Registries -# ------------------------------------------------------------------------------ - -CONDUCTOR_REGISTRY = {} +from bluemira.magnets.cable import ABCCable +from bluemira.magnets.utils import reciprocal_summation, summation -# ------------------------------------------------------------------------------ -# Strand Class -# ------------------------------------------------------------------------------ +if TYPE_CHECKING: + from matproplib import OperationalConditions + from matproplib.material import Material -class Conductor(metaclass=RegistrableMeta): +class Conductor: """ A generic conductor consisting of a cable surrounded by a jacket and an insulator. """ - _registry_ = CONDUCTOR_REGISTRY - _name_in_registry_ = "Conductor" - def __init__( self, cable: ABCCable, @@ -69,25 +51,35 @@ def __init__( mat_ins: insulator's material dx_jacket: - x-thickness of the jacket + x-thickness of the jacket [m]. dy_jacket: - y-tickness of the jacket + y-thickness of the jacket [m]. dx_ins: - x-thickness of the insulator + x-thickness of the insulator [m]. dy_ins: - y-tickness of the insulator + y-thickness of the insulator [m]. name: string identifier """ self.name = name - self._dx_jacket = dx_jacket + self.dx_jacket = dx_jacket self._dy_jacket = dy_jacket + self.dx_ins = dx_ins self._dy_ins = dy_ins - self._dx_ins = dx_ins self.mat_ins = mat_ins self.mat_jacket = mat_jacket self.cable = cable + @property + def dy_jacket(self): + """y-thickness of the jacket [m]""" + return self._dy_jacket + + @property + def dy_ins(self): + """y-thickness of the ins [m]""" + return self._dy_ins + @property def dx(self): """x-dimension of the conductor [m]""" @@ -98,42 +90,6 @@ def dy(self): """y-dimension of the conductor [m]""" return self.dy_ins * 2 + self.dy_jacket * 2 + self.cable.dy - @property - def dx_jacket(self): - """Thickness in the x-direction of the jacket [m]""" - return self._dx_jacket - - @dx_jacket.setter - def dx_jacket(self, value): - self._dx_jacket = value - - @property - def dy_jacket(self): - """Thickness in the y-direction of the jacket [m]""" - return self._dy_jacket - - @dy_jacket.setter - def dy_jacket(self, value): - self._dy_jacket = value - - @property - def dx_ins(self): - """Thickness in the x-direction of the insulator [m]""" - return self._dx_ins - - @dx_ins.setter - def dx_ins(self, value): - self._dx_ins = value - - @property - def dy_ins(self): - """Thickness in the y-direction of the jacket [m]""" - return self._dy_ins - - @dy_ins.setter - def dy_ins(self, value): - self._dy_ins = value - @property def area(self): """Area of the conductor [m^2]""" @@ -153,23 +109,21 @@ def area_ins(self): Returns ------- - float [m²] + : + area [m²] """ return self.area - self.area_jacket - self.cable.area - def to_dict(self) -> dict: + def to_dict(self) -> dict[str, Any]: """ - Serialize the conductor instance to a dictionary. + Serialise the conductor instance to a dictionary. Returns ------- - dict - Dictionary with serialized conductor data. + : + Dictionary with serialised conductor data. """ return { - "name_in_registry": getattr( - self, "_name_in_registry_", self.__class__.__name__ - ), "name": self.name, "cable": self.cable.to_dict(), "mat_jacket": self.mat_jacket, @@ -180,84 +134,21 @@ def to_dict(self) -> dict: "dy_ins": self.dy_ins, } - @classmethod - def from_dict( - cls, - conductor_dict: dict[str, Any], - name: str | None = None, - ) -> "Conductor": - """ - Deserialize a Conductor instance from a dictionary. - - Parameters - ---------- - cls : type - Class to instantiate (Conductor or subclass). - conductor_dict : dict - Dictionary containing serialized conductor data. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - Conductor - A fully reconstructed Conductor instance. - - Raises - ------ - ValueError - If the 'name_in_registry' field does not match the expected class - registration name, - or if the name already exists and unique_name is False. - """ - # Validate registration name - name_in_registry = conductor_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - # Deserialize cable - cable = create_cable_from_dict( - cable_dict=conductor_dict["cable"], - ) - - # Resolve jacket material - mat_jacket = conductor_dict["mat_jacket"] - - # Resolve insulation material - mat_ins = conductor_dict["mat_ins"] - - # Instantiate - return cls( - cable=cable, - mat_jacket=mat_jacket, - mat_ins=mat_ins, - dx_jacket=conductor_dict["dx_jacket"], - dy_jacket=conductor_dict["dy_jacket"], - dx_ins=conductor_dict["dx_ins"], - dy_ins=conductor_dict["dy_ins"], - name=name or conductor_dict.get("name"), - ) - - def erho(self, op_cond: OperationalConditions): + def erho(self, op_cond: OperationalConditions) -> float: """ Computes the conductor's equivalent resistivity considering the resistance of its strands in parallel. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float [Ohm m] + : + resistivity [Ohm m] Notes ----- @@ -267,23 +158,24 @@ def erho(self, op_cond: OperationalConditions): self.cable.erho(op_cond) / self.cable.area, self.mat_jacket.electrical_resistivity(op_cond) / self.area_jacket, ]) - res_tot = parall_r(resistances) + res_tot = reciprocal_summation(resistances) return res_tot * self.area - def Cp(self, op_cond: OperationalConditions): # noqa: N802 + def Cp(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Computes the conductor's equivalent specific heat considering the specific heats of its components in series. Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float [J/K/m] + : + Specific heat capacity [J/K/m] Notes ----- @@ -293,89 +185,72 @@ def Cp(self, op_cond: OperationalConditions): # noqa: N802 self.cable.Cp(op_cond) * self.cable.area, self.mat_jacket.specific_heat_capacity(op_cond) * self.area_jacket, ]) - return serie_r(weighted_specific_heat) / self.area - - def _mat_ins_y_modulus(self, op_cond: OperationalConditions): - return self.mat_ins.youngs_modulus(op_cond) - - def _mat_jacket_y_modulus(self, op_cond: OperationalConditions): - return self.mat_jacket.youngs_modulus(op_cond) + return summation(weighted_specific_heat) / self.area - def _Kx_topbot_ins(self, op_cond: OperationalConditions): # noqa: N802 + def _Kx_topbot_ins(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the top/bottom insulator in the x-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return self._mat_ins_y_modulus(op_cond) * self.cable.dy / self.dx_ins + return self.mat_ins.youngs_modulus(op_cond) * self.cable.dy / self.dx_ins - def _Kx_lat_ins(self, op_cond: OperationalConditions): # noqa: N802 + def _Kx_lat_ins(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the lateral insulator in the x-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return self._mat_ins_y_modulus(op_cond) * self.dy_ins / self.dx + return self.mat_ins.youngs_modulus(op_cond) * self.dy_ins / self.dx - def _Kx_lat_jacket(self, op_cond: OperationalConditions): # noqa: N802 + def _Kx_lat_jacket(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the lateral jacket in the x-direction. Returns ------- - float + : Axial stiffness [N/m] """ return ( - self._mat_jacket_y_modulus(op_cond) + self.mat_jacket.youngs_modulus(op_cond) * self.dy_jacket / (self.dx - 2 * self.dx_ins) ) - def _Kx_topbot_jacket(self, op_cond: OperationalConditions): # noqa: N802 + def _Kx_topbot_jacket(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the top/bottom jacket in the x-direction. Returns ------- - float - Axial stiffness [N/m] - """ - return self._mat_jacket_y_modulus(op_cond) * self.cable.dy / self.dx_jacket - - def _Kx_cable(self, op_cond: OperationalConditions): # noqa: N802 - """ - Equivalent stiffness of the cable in the x-direction. - - Returns - ------- - float + : Axial stiffness [N/m] """ - return self.cable.Kx(op_cond) + return self.mat_jacket.youngs_modulus(op_cond) * self.cable.dy / self.dx_jacket - def Kx(self, op_cond: OperationalConditions): # noqa: N802 + def Kx(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the conductor in the x-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return parall_k([ + return summation([ self._Kx_lat_ins(op_cond), self._Kx_lat_jacket(op_cond), - serie_k([ + reciprocal_summation([ self._Kx_topbot_ins(op_cond), self._Kx_topbot_jacket(op_cond), - self._Kx_cable(op_cond), + self.cable.Kx(op_cond), self._Kx_topbot_jacket(op_cond), self._Kx_topbot_ins(op_cond), ]), @@ -383,81 +258,70 @@ def Kx(self, op_cond: OperationalConditions): # noqa: N802 self._Kx_lat_ins(op_cond), ]) - def _Ky_topbot_ins(self, op_cond: OperationalConditions): # noqa: N802 + def _Ky_topbot_ins(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the top/bottom insulator in the y-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return self._mat_ins_y_modulus(op_cond) * self.cable.dx / self.dy_ins + return self.mat_ins.youngs_modulus(op_cond) * self.cable.dx / self.dy_ins - def _Ky_lat_ins(self, op_cond: OperationalConditions): # noqa: N802 + def _Ky_lat_ins(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the lateral insulator in the y-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return self._mat_ins_y_modulus(op_cond) * self.dx_ins / self.dy + return self.mat_ins.youngs_modulus(op_cond) * self.dx_ins / self.dy - def _Ky_lat_jacket(self, op_cond: OperationalConditions): # noqa: N802 + def _Ky_lat_jacket(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the lateral jacket in the y-direction. Returns ------- - float + : Axial stiffness [N/m] """ return ( - self._mat_jacket_y_modulus(op_cond) + self.mat_jacket.youngs_modulus(op_cond) * self.dx_jacket / (self.dy - 2 * self.dy_ins) ) - def _Ky_topbot_jacket(self, op_cond: OperationalConditions): # noqa: N802 + def _Ky_topbot_jacket(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the top/bottom jacket in the y-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return self._mat_jacket_y_modulus(op_cond) * self.cable.dx / self.dy_jacket + return self.mat_jacket.youngs_modulus(op_cond) * self.cable.dx / self.dy_jacket - def _Ky_cable(self, op_cond: OperationalConditions): # noqa: N802 - """ - Equivalent stiffness of the cable in the y-direction. - - Returns - ------- - float - Axial stiffness [N/m] - """ - return self.cable.Ky(op_cond) - - def Ky(self, op_cond: OperationalConditions): # noqa: N802 + def Ky(self, op_cond: OperationalConditions) -> float: # noqa: N802 """ Equivalent stiffness of the conductor in the y-direction. Returns ------- - float + : Axial stiffness [N/m] """ - return parall_k([ + return summation([ self._Ky_lat_ins(op_cond), self._Ky_lat_jacket(op_cond), - serie_k([ + reciprocal_summation([ self._Ky_topbot_ins(op_cond), self._Ky_topbot_jacket(op_cond), - self._Ky_cable(op_cond), + self.cable.Ky(op_cond), self._Ky_topbot_jacket(op_cond), self._Ky_topbot_ins(op_cond), ]), @@ -481,20 +345,21 @@ def _tresca_sigma_jacket( Parameters ---------- - pressure : + pressure: The pressure applied along the specified direction (Pa). - f_z : + f_z: The force applied in the z direction, perpendicular to the conductor cross-section (N). op_cond: OperationalConditions Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. - direction : + direction: The direction along which the pressure is applied ('x' or 'y'). Default is 'x'. Returns ------- + : The calculated Tresca stress in the jacket (Pa). Raises @@ -508,11 +373,11 @@ def _tresca_sigma_jacket( if direction == "x": saf_jacket = (self.cable.dx + 2 * self.dx_jacket) / (2 * self.dx_jacket) - K = parall_k([ # noqa: N806 + K = summation([ # noqa: N806 2 * self._Ky_lat_ins(op_cond), 2 * self._Ky_lat_jacket(op_cond), - serie_k([ - self._Ky_cable(op_cond), + reciprocal_summation([ + self.cable.Ky(op_cond), self._Ky_topbot_jacket(op_cond) / 2, ]), ]) @@ -522,210 +387,122 @@ def _tresca_sigma_jacket( else: saf_jacket = (self.cable.dy + 2 * self.dy_jacket) / (2 * self.dy_jacket) - K = parall_k([ # noqa: N806 + K = summation([ # noqa: N806 2 * self._Kx_lat_ins(op_cond), 2 * self._Kx_lat_jacket(op_cond), - serie_k([ - self._Kx_cable(op_cond), + reciprocal_summation([ + self.cable.Kx(op_cond), self._Kx_topbot_jacket(op_cond) / 2, ]), ]) X_jacket = 2 * self._Kx_lat_jacket(op_cond) / K # noqa: N806 - # tresca_stress = pressure * X_jacket * saf_jacket + f_z / self.area_jacket - + # tresca_stress return pressure * X_jacket * saf_jacket + f_z / self.area_jacket - def optimize_jacket_conductor( + def sigma_difference( self, + jacket_thickness: float, pressure: float, - f_z: float, + fz: float, op_cond: OperationalConditions, allowable_sigma: float, - bounds: np.ndarray | None = None, direction: str = "x", - ): + ) -> float: """ - Optimize the jacket dimension of a conductor based on allowable stress using - the Tresca criterion. + Objective function for optimising conductor jacket thickness based on the + Tresca yield criterion. + + This function computes the absolute difference between the calculated Tresca + stress in the jacket and the allowable stress. It is used as a fitness + function during scalar minimisation to determine the optimal jacket + thickness. Parameters ---------- - pressure : - The pressure applied along the specified direction (Pa). - f_z : - The force applied in the z direction, perpendicular to the conductor - cross-section (N). - op_cond: OperationalConditions + jacket_thickness: + Proposed thickness of the conductor jacket [m] in the direction + perpendicular to the applied pressure. + pressure: + Magnetic or mechanical pressure applied along the specified direction + [Pa]. + fz: + Axial or vertical force applied perpendicular to the cross-section [N]. + op_cond: Operational conditions including temperature, magnetic field, and strain - at which to calculate the material properties. - allowable_sigma : - The allowable stress (Pa) for the jacket material. - bounds : - Optional bounds for the jacket thickness optimization (default is None). - direction : - The direction along which the pressure is applied ('x' or 'y'). Default is - 'x'. + at which to calculate the material property. + allowable_sigma: + Maximum allowed stress for the jacket material [Pa]. + direction: + Direction of the applied pressure. Can be either 'x' (horizontal) or + 'y' (vertical). Default is 'x'. Returns ------- - The result of the optimization process containing information about the - optimal jacket thickness. + : + Absolute difference between the calculated Tresca stress and the + allowable stress [Pa]. Raises ------ ValueError - If the optimization process did not converge. + If the `direction` is not 'x' or 'y'. Notes ----- - This function uses the Tresca yield criterion to optimize the thickness of the - jacket surrounding the conductor. - This function directly update the conductor's jacket thickness along the x - direction to the optimal value. - """ - - def sigma_difference( - jacket_thickness: float, - pressure: float, - fz: float, - op_cond: OperationalConditions, - allowable_sigma: float, - direction: str = "x", - ): - """ - Objective function for optimizing conductor jacket thickness based on the - Tresca yield criterion. - - This function computes the absolute difference between the calculated Tresca - stress in the jacket and the allowable stress. It is used as a fitness - function during scalar minimization to determine the optimal jacket - thickness. - - Parameters - ---------- - jacket_thickness : float - Proposed thickness of the conductor jacket [m] in the direction - perpendicular to the applied pressure. - pressure : float - Magnetic or mechanical pressure applied along the specified direction - [Pa]. - fz : float - Axial or vertical force applied perpendicular to the cross-section [N]. - op_cond: OperationalConditions - Operational conditions including temperature, magnetic field, and strain - at which to calculate the material property. - allowable_sigma : float - Maximum allowed stress for the jacket material [Pa]. - direction : str, optional - Direction of the applied pressure. Can be either 'x' (horizontal) or - 'y' (vertical). Default is 'x'. - - Returns - ------- - float - Absolute difference between the calculated Tresca stress and the - allowable stress [Pa]. - - Raises - ------ - ValueError - If the `direction` is not 'x' or 'y'. - - Notes - ----- - - This function updates the conductor's internal jacket dimension ( - `dx_jacket` or `dy_jacket`) with the trial value `jacket_thickness`. - - It is intended for use with scalar optimization algorithms such as - `scipy.optimize.minimize_scalar`. - """ - if direction not in {"x", "y"}: - raise ValueError("Invalid direction: choose either 'x' or 'y'.") - - if direction == "x": - self.dx_jacket = jacket_thickness - else: - self.dy_jacket = jacket_thickness - - sigma_r = self._tresca_sigma_jacket(pressure, fz, op_cond, direction) - - # Normal difference - diff = abs(sigma_r - allowable_sigma) - - # Penalty if stress exceeds allowable - if sigma_r > allowable_sigma: - penalty = 1e6 + (sigma_r - allowable_sigma) * 1e6 - return diff + penalty - - return diff - - debug_msg = ["Method optimize_jacket_conductor:"] + - This function updates the conductor's internal jacket dimension ( + `dx_jacket` or `dy_jacket`) with the trial value `jacket_thickness`. + - It is intended for use with scalar optimisation algorithms such as + `scipy.optimize.minimize_scalar`. + """ + if direction not in {"x", "y"}: + raise ValueError("Invalid direction: choose either 'x' or 'y'.") if direction == "x": - debug_msg.append(f"Previous dx_jacket: {self.dx_jacket}") + self.dx_jacket = jacket_thickness else: - debug_msg.append(f"Previous dy_jacket: {self.dy_jacket}") + self.dy_jacket = jacket_thickness - method = "bounded" if bounds is not None else None + sigma_r = self._tresca_sigma_jacket(pressure, fz, op_cond, direction) - if method == "bounded": - debug_msg.append(f"bounds: {bounds}") + # Normal difference + diff = abs(sigma_r - allowable_sigma) - result = minimize_scalar( - fun=sigma_difference, - args=(pressure, f_z, op_cond, allowable_sigma), - bounds=bounds, - method=method, - options={"xatol": 1e-4}, - ) + # Penalty if stress exceeds allowable + if sigma_r > allowable_sigma: + penalty = 1e6 + (sigma_r - allowable_sigma) * 1e6 + return diff + penalty - if not result.success: - raise ValueError("Optimization of the jacket conductor did not converge.") - if direction == "x": - self.dx_jacket = result.x - debug_msg.append(f"Optimal dx_jacket: {self.dx_jacket}") - else: - self.dy_jacket = result.x - debug_msg.append(f"Optimal dy_jacket: {self.dy_jacket}") - debug_msg.append( - f"Averaged sigma in the {direction}-direction: " - f"{self._tresca_sigma_jacket(pressure, f_z, op_cond) / 1e6} MPa\n" - f"Allowable stress in the {direction}-direction: {allowable_sigma / 1e6} " - f"MPa." - ) - debug_msg = "\n".join(debug_msg) - bluemira_debug(debug_msg) - - return result + return diff def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): """ Plot a schematic cross-section of the conductor, including cable, jacket, and insulator layers. - This method visualizes the hierarchical geometry of the conductor centered + This method visualises the hierarchical geometry of the conductor centered at a given position. The jacket and insulator are drawn as rectangles, while the internal cable uses its own plotting method. Parameters ---------- - xc : float, optional + xc: X-coordinate of the conductor center in the reference coordinate system. Default is 0. - yc : float, optional + yc: Y-coordinate of the conductor center in the reference coordinate system. Default is 0. - show : bool, optional + show: If True, the figure is rendered immediately using `plt.show()`. Default is False. - ax : matplotlib.axes.Axes or None, optional + ax: Axis on which to render the plot. If None, a new figure and axis will be created internally. Returns ------- - ax : matplotlib.axes.Axes + ax: The axis containing the rendered plot. Notes @@ -770,13 +547,13 @@ def plot(self, xc: float = 0, yc: float = 0, *, show: bool = False, ax=None): return ax - def __str__(self): + def __str__(self) -> str: """ Generate a human-readable string representation of the conductor. Returns ------- - str + : A multi-line summary of the conductor's key dimensions and its nested cable description. This includes: - Total x and y dimensions, @@ -803,8 +580,6 @@ class SymmetricConductor(Conductor): mantain a constant thickness (i.e. dy_jacket = dx_jacket and dy_ins = dx_ins). """ - _name_in_registry_ = "SymmetricConductor" - def __init__( self, cable: ABCCable, @@ -827,9 +602,9 @@ def __init__( mat_ins: insulator's material dx_jacket: - x(y)-thickness of the jacket + x-thickness of the jacket [m]. dx_ins: - x(y)-thickness of the insulator + x-thickness of the insulator [m]. name: string identifier @@ -848,13 +623,13 @@ def __init__( ) @property - def dy_jacket(self): + def dy_jacket(self) -> float: """ y-thickness of the jacket [m]. Returns ------- - float + : Notes ----- @@ -863,13 +638,13 @@ def dy_jacket(self): return self.dx_jacket @property - def dy_ins(self): + def dy_ins(self) -> float: """ y-thickness of the insulator [m]. Returns ------- - float + : Notes ----- @@ -879,123 +654,18 @@ def dy_ins(self): def to_dict(self) -> dict: """ - Serialize the symmetric conductor instance to a dictionary. + Serialise the symmetric conductor instance to a dictionary. Returns ------- - dict - Dictionary with serialized symmetric conductor data. + : + Dictionary with serialised symmetric conductor data. """ return { - "name_in_registry": getattr( - self, "_name_in_registry_", self.__class__.__name__ - ), "name": self.name, "cable": self.cable.to_dict(), - "mat_jacket": self.mat_jacket, - "mat_ins": self.mat_ins, + "mat_jacket": self.mat_jacket.name, + "mat_ins": self.mat_ins.name, "dx_jacket": self.dx_jacket, "dx_ins": self.dx_ins, } - - @classmethod - def from_dict( - cls, - conductor_dict: dict[str, Any], - name: str | None = None, - ) -> "SymmetricConductor": - """ - Deserialize a SymmetricConductor instance from a dictionary. - - Parameters - ---------- - cls : type - Class to instantiate (SymmetricConductor). - conductor_dict : dict - Dictionary containing serialized conductor data. - name : str, optional - Name for the new instance. - - Returns - ------- - SymmetricConductor - A fully reconstructed SymmetricConductor instance. - - Raises - ------ - ValueError - If the 'name_in_registry' does not match the expected registration name. - """ - # Validate registration name - name_in_registry = conductor_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - # Deserialize cable - cable = create_cable_from_dict( - cable_dict=conductor_dict["cable"], - ) - - # Resolve jacket material - mat_jacket = conductor_dict["mat_jacket"] - - # Resolve insulation material - mat_ins = conductor_dict["mat_ins"] - - # Instantiate - return cls( - cable=cable, - mat_jacket=mat_jacket, - mat_ins=mat_ins, - dx_jacket=conductor_dict["dx_jacket"], - dx_ins=conductor_dict["dx_ins"], - name=name or conductor_dict.get("name"), - ) - - -def create_conductor_from_dict( - conductor_dict: dict, - name: str | None = None, -) -> "Conductor": - """ - Factory function to create a Conductor (or subclass) from a serialized dictionary. - - Parameters - ---------- - conductor_dict : dict - Serialized conductor dictionary, must include 'name_in_registry' field. - name : str, optional - Name to assign to the created conductor. If None, uses the name in the - dictionary. - - Returns - ------- - Conductor - A fully instantiated Conductor (or subclass) object. - - Raises - ------ - ValueError - If no class is registered with the given name_in_registry. - """ - name_in_registry = conductor_dict.get("name_in_registry") - if name_in_registry is None: - raise ValueError("Conductor dictionary must include 'name_in_registry' field.") - - conductor_cls = CONDUCTOR_REGISTRY.get(name_in_registry) - if conductor_cls is None: - available = list(CONDUCTOR_REGISTRY.keys()) - raise ValueError( - f"No registered conductor class with name_in_registry '{name_in_registry}'. " - f"Available: {available}" - ) - - return conductor_cls.from_dict( - name=name, - conductor_dict=conductor_dict, - ) diff --git a/bluemira/magnets/fatigue.py b/bluemira/magnets/fatigue.py index 3962f31031..c4bbd4b528 100644 --- a/bluemira/magnets/fatigue.py +++ b/bluemira/magnets/fatigue.py @@ -8,8 +8,11 @@ Paris Law fatigue model with FE-inspired analytical crack propagation """ +from __future__ import annotations + import abc from dataclasses import dataclass +from typing import Final import numpy as np @@ -87,7 +90,7 @@ def _boundary_correction_factor( Returns ------- - float + : Boundary correction factor F. """ return (m1 + m2 * a_d_t**2 + m3 * a_d_t**4) * g * f_phi * f_w @@ -99,7 +102,7 @@ def _bending_correction_factor(h1: float, h2: float, p: float, phi: float) -> fl Returns ------- - float + : Bending correction factor. """ return h1 + (h2 - h1) * np.sin(phi) ** p @@ -111,7 +114,7 @@ def _ellipse_shape_factor(ratio: float) -> float: Returns ------- - float + : Shape factor Q. """ return 1.0 + 1.464 * ratio**1.65 @@ -123,7 +126,7 @@ def _angular_location_correction(a: float, c: float, phi: float) -> float: Returns ------- - float + : Angular correction factor f_phi. """ if a <= c: @@ -137,7 +140,7 @@ def _finite_width_correction(a_d_t: float, c: float, w: float) -> float: Returns ------- - float + : Finite width correction factor. """ return 1.0 / np.sqrt(np.cos(np.sqrt(a_d_t) * np.pi * c / (2 * w))) # (11) @@ -155,20 +158,18 @@ class Crack(abc.ABC): Crack width along the plate length direction """ - alpha = None - - def __init__(self, depth: float, width: float): - self.depth = depth # a - self.width = width # c + def __init__(self, depth, width): + self.depth = depth + self.width = width @classmethod - def from_area(cls, area: float, aspect_ratio: float): + def from_area(cls, area: float, aspect_ratio: float) -> Crack: """ Instatiate a crack from an area and aspect ratio Returns ------- - Crack + : New instance of the crack geometry. """ depth = np.sqrt(area / (cls.alpha * np.pi * aspect_ratio)) @@ -182,11 +183,15 @@ def area(self) -> float: Returns ------- - float + : Area [m²]. """ return self.alpha * np.pi * self.depth * self.width + @property + @abc.abstractmethod + def alpha(self) -> float: ... + @abc.abstractmethod def stress_intensity_factor( self, @@ -215,7 +220,7 @@ class QuarterEllipticalCornerCrack(Crack): Crack width along the plate length direction """ - alpha = 0.25 + alpha: Final[float] = 0.25 def stress_intensity_factor( # noqa: PLR6301 self, @@ -249,7 +254,8 @@ def stress_intensity_factor( # noqa: PLR6301 Returns ------- - Stress intensity factor + : + Stress intensity factor Notes ----- @@ -315,7 +321,7 @@ class SemiEllipticalSurfaceCrack(Crack): Crack width along the plate length direction """ - alpha = 0.5 + alpha: Final[float] = 0.5 def stress_intensity_factor( # noqa: PLR6301 self, @@ -349,7 +355,8 @@ def stress_intensity_factor( # noqa: PLR6301 Returns ------- - Stress intensity factor + : + Stress intensity factor Notes ----- @@ -406,7 +413,7 @@ class EllipticalEmbeddedCrack(Crack): Crack width along the plate length direction """ - alpha = 1.0 + alpha: Final[float] = 1.0 def stress_intensity_factor( # noqa: PLR6301 self, @@ -440,7 +447,8 @@ def stress_intensity_factor( # noqa: PLR6301 Returns ------- - Stress intensity factor + : + Stress intensity factor Notes ----- @@ -495,7 +503,8 @@ def calculate_n_pulses( Returns ------- - Number of plasma pulses + : + Number of plasma pulses Notes ----- diff --git a/bluemira/magnets/init_magnets_registry.py b/bluemira/magnets/init_magnets_registry.py deleted file mode 100644 index 5f47977230..0000000000 --- a/bluemira/magnets/init_magnets_registry.py +++ /dev/null @@ -1,68 +0,0 @@ -# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza -# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh -# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short -# -# SPDX-License-Identifier: LGPL-2.1-or-later - -""" -Initialization functions to register magnet classes. - -These functions import necessary modules to trigger class registration -(via metaclasses) without polluting the importing namespace. -""" - - -def register_strands(): - """ - Import and register all known Strand classes. - - This triggers their metaclass registration into the STRAND_REGISTRY. - Importing here avoids polluting the top-level namespace. - - Classes registered - ------------------- - - Strand - - SuperconductingStrand - """ - - -def register_cables(): - """ - Import and register all known Cable classes. - - This triggers their metaclass registration into the CABLE_REGISTRY. - Importing here avoids polluting the top-level namespace. - - Classes registered - ------------------- - - Cable - - Specialized cable types (e.g., TwistedCables) - """ - - -def register_conductors(): - """ - Import and register all known Conductor classes. - - This triggers their metaclass registration into the CONDUCTOR_REGISTRY. - Importing here avoids polluting the top-level namespace. - - Classes registered - ------------------- - - Conductor - - Specialized conductors - """ - - -def register_all_magnets(): - """ - Import and register all known magnet-related classes. - - Calls `register_strands()`, `register_cables()`, and `register_conductors()` - to fully populate all internal registries. - - Use this function at initialization if you want all classes to be available. - """ - register_strands() - register_cables() - register_conductors() diff --git a/bluemira/magnets/registry.py b/bluemira/magnets/registry.py deleted file mode 100644 index 8159ab3acd..0000000000 --- a/bluemira/magnets/registry.py +++ /dev/null @@ -1,445 +0,0 @@ -# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza -# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh -# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short -# -# SPDX-License-Identifier: LGPL-2.1-or-later -""" -Generic class and instance registration utilities. - -This module provides: -- RegistrableMeta: A metaclass that automatically registers classes into a specified -registry. -- InstanceRegistrable: A mixin that automatically registers instances into a specified -global cache. - -Intended for use in frameworks where automatic discovery of classes -and instances is required, such as for strands, cables, conductors, or other physical -models. - -Usage ------ -Classes intended to be registered must: -- Define a class-level `_registry_` dictionary (for class registration). -- Optionally set a `_name_in_registry_` string (custom name for registration). - -Instances intended to be globally tracked must: -- Inherit from InstanceRegistrable. -- Provide a unique `name` attribute at creation. -""" - -from abc import ABCMeta -from typing import ClassVar - -from bluemira.base.look_and_feel import bluemira_debug - - -# ------------------------------------------------------------------------------ -# RegistrableMeta -# ------------------------------------------------------------------------------ -class RegistrableMeta(ABCMeta): - """ - Metaclass for automatic class registration into a registry. - - Enforces that: - - '_name_in_registry_' must be explicitly defined in every class body (no - inheritance allowed). - - '_registry_' can be inherited if not redefined. - """ - - def __new__(mcs, name, bases, namespace): - """ - Create and register a new class instance using the RegistrableMeta metaclass. - - This method: - - Automatically registers concrete (non-abstract) classes into a specified - registry. - - Enforces that concrete classes explicitly declare a '_name_in_registry_'. - - Allows '_registry_' to be inherited from base classes if not redefined. - - Parameters - ---------- - mcs : type - The metaclass (usually RegistrableMeta itself). - name : str - The name of the class being created. - bases : tuple of type - The base classes of the class being created. - namespace : dict - The attribute dictionary of the class. - - Returns - ------- - type - The newly created class. - - Raises - ------ - TypeError - If a concrete class does not define a '_name_in_registry_'. - If no '_registry_' can be found (either defined or inherited). - ValueError - If a duplicate '_name_in_registry_' is detected within the registry. - - Notes - ----- - Abstract base classes (ABCs) are exempted from registration requirements. - - Registration process: - - If the class is abstract, skip registration. - - Otherwise: - - Check for existence of '_registry_' (allow inheritance). - - Require explicit '_name_in_registry_' (must be defined in the class body). - - Insert the class into the registry under the specified name. - """ - bluemira_debug(f"Registering {name}...") # Debug print - - cls = super().__new__(mcs, name, bases, namespace) - - is_abstract = bool(getattr(cls, "__abstractmethods__", False)) - - if not is_abstract: - # Only enforce _name_in_registry_ and _registry_ for concrete classes - # _registry_ can be inherited - registry = getattr(cls, "_registry_", None) - - # _name_in_registry_ must be explicit in the class body - register_name = namespace.get("_name_in_registry_", None) - - # Checks - if registry is None: - raise TypeError( - f"Class {name} must define or inherit a '_registry_' for " - f"registration." - ) - - if register_name is None: - raise TypeError( - f"Class {name} must explicitly define a '_name_in_registry_' for " - f"registration." - ) - - # Registration - if register_name: - if register_name in registry: - raise ValueError( - f"Duplicate registration for class '{register_name}'." - ) - registry[register_name] = cls - cls._name_in_registry_ = register_name # Optional: for introspection - - return cls - - @classmethod - def unregister(cls): - """ - Unregister the class from its associated registry. - - This method removes the class from the `_registry_` dictionary under its - '_name_in_registry_' key. It is safe to call at runtime to dynamically - de-register classes, such as when reloading modules or cleaning up. - - Raises - ------ - AttributeError - If the class does not have a '_registry_' or '_name_in_registry_' attribute. - KeyError - If the class is not found in the registry (already unregistered or never - registered). - - Notes - ----- - - Only registered (non-abstract) classes have the 'unregister' method. - - Abstract base classes (ABCs) are skipped and do not perform registration. - """ - registry = getattr(cls, "_registry_", None) - name_in_registry = getattr(cls, "_name_in_registry_", None) - - if registry is None or name_in_registry is None: - raise AttributeError( - "Cannot unregister: missing '_registry_' or '_name_in_registry_'." - ) - - try: - del registry[name_in_registry] - except KeyError: - raise KeyError( - f"Class '{name_in_registry}' not found in the registry; it may have " - f"already been unregistered." - ) from None - - @classmethod - def get_registered_class(cls, name: str): - """ - Retrieve a registered class by name. - - Parameters - ---------- - name : str - Name of the registered class. - - Returns - ------- - type or None - The registered class, or None if not found. - - Raises - ------ - AttributeError - If the class does not define or inherit a '_registry_' attribute. - """ - registry = getattr(cls, "_registry_", None) - if registry is None: - raise AttributeError( - f"Class {cls.__name__} must define or inherit a '_registry_' attribute." - ) - return registry.get(name) - - @classmethod - def list_registered_classes(cls) -> list[str]: - """ - List names of all registered classes. - - Returns - ------- - list of str - List of names of registered classes. - - Raises - ------ - AttributeError - If the class does not define or inherit a '_registry_' attribute. - """ - registry = getattr(cls, "_registry_", None) - if registry is None: - raise AttributeError( - f"Class {cls.__name__} must define or inherit a '_registry_' attribute." - ) - return list(registry.keys()) - - @classmethod - def clear_registered_classes(cls): - """ - Clear all registered classes from the registry. - - This method removes all entries from the `_registry_` dictionary, - effectively unregistering all previously registered classes. - - Raises - ------ - AttributeError - If the class does not define or inherit a '_registry_' attribute. - """ - registry = getattr(cls, "_registry_", None) - if registry is None: - raise AttributeError( - f"Class {cls.__name__} must define or inherit a '_registry_' attribute." - ) - registry.clear() - - -# ------------------------------------------------------------------------------ -# InstanceRegistrable -# ------------------------------------------------------------------------------ - - -class InstanceRegistrable: - """ - Mixin class to automatically register instances into a global instance cache. - - This class provides: - - Automatic instance registration into a global cache. - - Optional control over registration (register or not). - - Optional automatic generation of unique names to avoid conflicts. - - Attributes - ---------- - _global_instance_cache_ : dict - Class-level cache shared among all instances for lookup by name. - _do_not_register : bool - If True, the instance will not be registered. - _unique : bool - If True, automatically generate a unique name if the desired name already exists. - """ - - _global_instance_cache_: ClassVar[dict] = {} - - def __init__(self, name: str, *, unique_name: bool = False): - """ - Initialize an instance and optionally register it. - - Parameters - ---------- - name : str - Desired name of the instance. - unique_name : bool, optional - If True, generate a unique name if the given name already exists. - If False (strict mode), raise a ValueError on duplicate names. - """ - self._unique_name = None - self.unique_name = unique_name - - # Setting the name will trigger registration (unless do_registration is False) - self._name = None - self.name = name - - @property - def unique_name(self) -> bool: - """ - Flag indicating whether to automatically generate a unique name on conflict. - - Returns - ------- - bool - True if automatic unique name generation is enabled (the passed name is - neglected) - False if strict name checking is enforced. - """ - return self._unique_name - - @unique_name.setter - def unique_name(self, value: bool): - """ - Set whether automatic unique name generation should be enabled. - - Parameters - ---------- - value : bool - If True, automatically generate a unique name if the desired name - is already registered. - If False, raise a ValueError if the name already exists. - """ - self._unique_name = value - - @property - def name(self) -> str: - """Return the instance name.""" - return self._name - - @name.setter - def name(self, value: str): - """ - Set the instance name and (re)register it according to registration rules. - - Behavior - -------- - - If `_do_not_register` is True, just assign the name without caching. - - If `unique_name` is True and the name already exists, automatically generate - a unique name. - - If `unique_name` is False and the name already exists, raise ValueError. - - Parameters - ---------- - value : str - Desired instance name. - - Raises - ------ - ValueError - If `unique_name` is False and the name is already registered. - """ - if hasattr(self, "_name") and self._name is not None: - self._unregister_self() - - if value is None: - self._name = None - return - - if value in self._global_instance_cache_: - if self.unique_name: - value = self.generate_unique_name(value) - else: - raise ValueError(f"Instance with name '{value}' already registered.") - - self._name = value - self._register_self() - - def _register_self(self): - """ - Register this instance into the global instance cache. - - Raises - ------ - AttributeError - If the instance does not have a 'name' attribute. - ValueError - If an instance with the same name already exists and unique is False. - """ - if getattr(self, "_do_not_register", False): - return # Skip registration if explicitly disabled - - if not hasattr(self, "name") or self.name is None: - raise AttributeError("Instance must have a 'name' attribute to register.") - - if self.name in self._global_instance_cache_: - if self.unique_name: - self.name = self.generate_unique_name(self.name) - else: - raise ValueError(f"Instance with name '{self.name}' already registered.") - - self._global_instance_cache_[self.name] = self - - def _unregister_self(self): - """ - Unregister this instance from the global instance cache. - """ - if hasattr(self, "name") and self.name in self._global_instance_cache_: - del self._global_instance_cache_[self.name] - - @classmethod - def get_registered_instance(cls, name: str): - """ - Retrieve a registered instance by name. - - Parameters - ---------- - name : str - Name of the registered instance. - - Returns - ------- - InstanceRegistrable or None - The registered instance, or None if not found. - """ - return cls._global_instance_cache_.get(name) - - @classmethod - def list_registered_instances(cls) -> list[str]: - """ - List names of all registered instances. - - Returns - ------- - list of str - List of names of registered instances. - """ - return list(cls._global_instance_cache_.keys()) - - @classmethod - def clear_registered_instances(cls): - """ - Clear all registered instances from the global cache. - """ - cls._global_instance_cache_.clear() - - @classmethod - def generate_unique_name(cls, base_name: str) -> str: - """ - Generate a unique name by appending a numeric suffix if necessary. - - Parameters - ---------- - base_name : str - Desired base name. - - Returns - ------- - str - Unique name guaranteed not to conflict with existing instances. - """ - if base_name not in cls._global_instance_cache_: - return base_name - - i = 1 - while f"{base_name}_{i}" in cls._global_instance_cache_: - i += 1 - return f"{base_name}_{i}" diff --git a/bluemira/magnets/strand.py b/bluemira/magnets/strand.py index 8086f665aa..4b99bc3773 100644 --- a/bluemira/magnets/strand.py +++ b/bluemira/magnets/strand.py @@ -12,6 +12,8 @@ - Automatic class and instance registration mechanisms """ +from __future__ import annotations + from typing import Any import matplotlib.pyplot as plt @@ -24,61 +26,44 @@ from bluemira.display.plotter import PlotOptions from bluemira.geometry.face import BluemiraFace from bluemira.geometry.tools import make_circle -from bluemira.magnets.registry import RegistrableMeta - -# ------------------------------------------------------------------------------ -# Global Registries -# ------------------------------------------------------------------------------ - -STRAND_REGISTRY = {} -# ------------------------------------------------------------------------------ -# Strand Class -# ------------------------------------------------------------------------------ - -class Strand(metaclass=RegistrableMeta): +class Strand: """ - Represents a strand with a circular cross-section, composed of a homogenized + Represents a strand with a circular cross-section, composed of a homogenised mixture of materials. This class automatically registers itself and its instances. """ - _registry_ = STRAND_REGISTRY - _name_in_registry_ = "Strand" - def __init__( self, materials: list[MaterialFraction], - d_strand: float = 0.82e-3, - temperature: float | None = None, + d_strand: float, + operating_temperature: float, name: str | None = "Strand", ): """ - Initialize a Strand instance. + Initialise a Strand instance. Parameters ---------- - materials : list of MaterialFraction + materials: Materials composing the strand with their fractions. - d_strand : float, optional - Strand diameter in meters (default 0.82e-3). - temperature : float, optional + d_strand: + Strand diameter [m]. + operating_temperature: float Operating temperature [K]. - name : str or None, optional + + name: Name of the strand. Defaults to "Strand". """ - self._d_strand = None - self._shape = None - self._materials = None - self._temperature = None - self.d_strand = d_strand + self.operating_temperature = operating_temperature self.materials = materials - self.name = name - self.temperature = temperature + self.name = name + self._shape = None # Create homogenised material self._homogenised_material = mixture( name=name, @@ -87,25 +72,25 @@ def __init__( ) @property - def materials(self) -> list: + def materials(self) -> list[MaterialFraction]: """ List of MaterialFraction materials composing the strand. Returns ------- - list of MaterialFraction + : Materials and their fractions. """ return self._materials @materials.setter - def materials(self, new_materials: list): + def materials(self, new_materials: list[MaterialFraction]): """ Set a new list of materials for the strand. Parameters ---------- - new_materials : list of MaterialFraction + new_materials: New materials to set. Raises @@ -127,84 +112,6 @@ def materials(self, new_materials: list): self._materials = new_materials - @property - def temperature(self) -> float | None: - """ - Operating temperature of the strand. - - Returns - ------- - float or None - Temperature in Kelvin. - """ - return self._temperature - - @temperature.setter - def temperature(self, value: float | None): - """ - Set a new operating temperature for the strand. - - Parameters - ---------- - value : float or None - New operating temperature in Kelvin. - - Raises - ------ - ValueError - If temperature is negative. - TypeError - If temperature is not a float or None. - """ - if value is not None: - if not isinstance(value, (float, int)): - raise TypeError( - f"temperature must be a float or int, got {type(value).__name__}." - ) - - if value < 0: - raise ValueError("Temperature cannot be negative.") - - self._temperature = float(value) if value is not None else None - - @property - def d_strand(self) -> float: - """ - Diameter of the strand. - - Returns - ------- - Parameter - Diameter [m]. - """ - return self._d_strand - - @d_strand.setter - def d_strand(self, d: float): - """ - Set the strand diameter and reset shape if changed. - - Parameters - ---------- - d : float or Parameter - New strand diameter. - - Raises - ------ - ValueError - If diameter is non-positive. - TypeError - If diameter is not a float number. - """ - if not isinstance(d, (float, int)): - raise TypeError(f"d_strand must be a float, got {type(d).__name__}") - if d <= 0: - raise ValueError("d_strand must be positive.") - - if self.d_strand is None or d != self.d_strand: - self._d_strand = float(d) - self._shape = None - @property def area(self) -> float: """ @@ -212,7 +119,7 @@ def area(self) -> float: Returns ------- - float + : Area [m²]. """ return np.pi * (self.d_strand**2) / 4 @@ -224,7 +131,7 @@ def shape(self) -> BluemiraFace: Returns ------- - BluemiraFace + : Circular face of the strand. """ if self._shape is None: @@ -237,13 +144,13 @@ def E(self, op_cond: OperationalConditions) -> float: # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Young's modulus [Pa]. """ return self._homogenised_material.youngs_modulus(op_cond) @@ -254,13 +161,13 @@ def rho(self, op_cond: OperationalConditions) -> float: Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Density [kg/m³]. """ return self._homogenised_material.density(op_cond) @@ -271,13 +178,13 @@ def erho(self, op_cond: OperationalConditions) -> float: Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Electrical resistivity [Ohm·m]. """ # Treat parallel calculation for resistivity @@ -295,13 +202,13 @@ def Cp(self, op_cond: OperationalConditions) -> float: # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Specific heat [J/kg/K]. """ # Treat volume/specific heat capacity calculation @@ -320,17 +227,19 @@ def Cp(self, op_cond: OperationalConditions) -> float: # noqa: N802 ) return self._homogenised_material.specific_heat_capacity(op_cond) - def plot(self, ax=None, *, show: bool = True, **kwargs): + def plot( + self, ax: plt.Axes | None = None, *, show: bool = True, **kwargs + ) -> plt.Axes: """ Plot a 2D cross-section of the strand. Parameters ---------- - ax : matplotlib.axes.Axes, optional + ax: Axis to plot on. - show : bool, optional + show: Whether to show the plot immediately. - kwargs : dict + kwargs: Additional arguments passed to the plot function. Returns @@ -350,7 +259,7 @@ def __str__(self) -> str: Returns ------- - str + : Description of the strand. """ return ( @@ -360,22 +269,19 @@ def __str__(self) -> str: f"shape = {self.shape}\n" ) - def to_dict(self) -> dict: + def to_dict(self) -> dict[str, Any]: """ - Serialize the strand instance to a dictionary. + Serialise the strand instance to a dictionary. Returns ------- - dict - Dictionary with serialized strand data. + : + Dictionary with serialised strand data. """ return { - "name_in_registry": getattr( - self, "_name_in_registry_", self.__class__.__name__ - ), "name": self.name, "d_strand": self.d_strand, - "temperature": self.temperature, + "temperature": self.operating_temperature, "materials": [ { "material": m.material, @@ -385,76 +291,6 @@ def to_dict(self) -> dict: ], } - @classmethod - def from_dict( - cls, - strand_dict: dict[str, Any], - name: str | None = None, - ) -> "Strand": - """ - Deserialize a Strand instance from a dictionary. - - Parameters - ---------- - cls : type - Class to instantiate (Strand or subclass). - strand_dict : dict - Dictionary containing serialized strand data. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - Strand - A new instantiated Strand object. - - Raises - ------ - TypeError - If the materials in the dictionary are not valid MaterialFraction instances. - ValueError - If the name_in_registry in the dictionary does not match the expected - class registration name. - """ - # Validate registration name - name_in_registry = strand_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. " - f"Expected '{expected_name_in_registry}'." - ) - - # Deserialize materials - material_mix = [] - for m in strand_dict["materials"]: - material_data = m["material"] - if isinstance(material_data, str): - raise TypeError( - "Material data must be a Material instance, not a string - " - "TEMPORARY." - ) - material_obj = material_data - - material_mix.append( - MaterialFraction(material=material_obj, fraction=m["fraction"]) - ) - - return cls( - materials=material_mix, - temperature=strand_dict.get("temperature"), - d_strand=strand_dict.get("d_strand"), - name=name or strand_dict.get("name"), - ) - - -# ------------------------------------------------------------------------------ -# SuperconductingStrand Class -# ------------------------------------------------------------------------------ - class SuperconductingStrand(Strand): """ @@ -466,34 +302,32 @@ class SuperconductingStrand(Strand): Automatically registered using the RegistrableMeta metaclass. """ - _name_in_registry_ = "SuperconductingStrand" - def __init__( self, materials: list[MaterialFraction], - d_strand: float = 0.82e-3, - temperature: float | None = None, + d_strand: float, + operating_temperature: float, name: str | None = "SuperconductingStrand", ): """ - Initialize a superconducting strand. + Initialise a superconducting strand. Parameters ---------- - materials : list of MaterialFraction + materials: Materials composing the strand with their fractions. One material must be a supercoductor. - d_strand : float, optional - Strand diameter in meters (default 0.82e-3). - temperature : float, optional + d_strand: + Strand diameter [m]. + operating_temperature: float Operating temperature [K]. - name : str or None, optional + name: Name of the strand. Defaults to "Strand". """ super().__init__( materials=materials, d_strand=d_strand, - temperature=temperature, + operating_temperature=operating_temperature, name=name, ) self._sc = self._check_materials() @@ -504,7 +338,7 @@ def _check_materials(self) -> MaterialFraction: Returns ------- - MaterialFraction + : The identified superconducting material. Raises @@ -539,7 +373,7 @@ def sc_area(self) -> float: Returns ------- - float + : Superconducting area [m²]. """ return self.area * self._sc.fraction @@ -550,13 +384,13 @@ def Jc(self, op_cond: OperationalConditions) -> float: # noqa:N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Critical current density [A/m²]. """ if op_cond.strain is None: @@ -569,13 +403,13 @@ def Ic(self, op_cond: OperationalConditions) -> float: # noqa:N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Critical current [A]. """ return self.Jc(op_cond) * self.sc_area @@ -584,30 +418,30 @@ def plot_Ic_B( # noqa:N802 self, B: np.ndarray, temperature: float, - ax=None, + ax: plt.Axes | None = None, *, show: bool = True, **kwargs, - ): + ) -> plt.Axes: """ Plot critical current Ic as a function of magnetic field B. Parameters ---------- - B : np.ndarray + B: Array of magnetic field values [T]. - temperature : float + temperature: Operating temperature [K]. - ax : matplotlib.axes.Axes, optional + ax: Axis to plot on. If None, a new figure is created. - show : bool, optional + show: Whether to immediately show the plot. - kwargs : dict + kwargs: Additional arguments passed to Ic calculation. Returns ------- - matplotlib.axes.Axes + : Axis with the plotted Ic vs B curve. """ if ax is None: @@ -635,48 +469,3 @@ def plot_Ic_B( # noqa:N802 plt.show() return ax - - -# ------------------------------------------------------------------------------ -# Supporting functions -# ------------------------------------------------------------------------------ -def create_strand_from_dict( - strand_dict: dict[str, Any], - name: str | None = None, -): - """ - Factory function to create a Strand or its subclass from a serialized dictionary. - - Parameters - ---------- - strand_dict : dict - Dictionary with serialized strand data. Must include a 'name_in_registry' field - corresponding to a registered class. - name : str, optional - If given, overrides the name from the dictionary. - - Returns - ------- - Strand - An instance of the appropriate Strand subclass. - - Raises - ------ - ValueError - If 'name_in_registry' is missing from the dictionary. - If no matching registered class is found. - """ - name_in_registry = strand_dict.get("name_in_registry") - if name_in_registry is None: - raise ValueError( - "Serialized strand dictionary must contain a 'name_in_registry' field." - ) - - cls = STRAND_REGISTRY.get(name_in_registry) - if cls is None: - raise ValueError( - f"No registered strand class with registration name '{name_in_registry}'. " - "Available classes are: " + ", ".join(STRAND_REGISTRY.keys()) - ) - - return cls.from_dict(name=name, strand_dict=strand_dict) diff --git a/bluemira/magnets/tfcoil_designer.py b/bluemira/magnets/tfcoil_designer.py new file mode 100644 index 0000000000..79a358f2d5 --- /dev/null +++ b/bluemira/magnets/tfcoil_designer.py @@ -0,0 +1,1060 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later +"""Designer for TF Coil XY cross section.""" + +from __future__ import annotations + +from collections.abc import Callable, Sequence +from dataclasses import dataclass +from typing import TYPE_CHECKING, Any + +import matplotlib.pyplot as plt +import numpy as np +import numpy.typing as npt +from matproplib import OperationalConditions +from matproplib.material import MaterialFraction +from scipy.optimize import minimize_scalar + +from bluemira.base.constants import MU_0, MU_0_2PI, MU_0_4PI +from bluemira.base.designer import Designer +from bluemira.base.look_and_feel import ( + bluemira_debug, + bluemira_error, + bluemira_print, + bluemira_warn, +) +from bluemira.base.parameter_frame import Parameter, ParameterFrame +from bluemira.magnets.cable import ABCCable, RectangularCable +from bluemira.magnets.conductor import Conductor, SymmetricConductor +from bluemira.magnets.utils import delayed_exp_func +from bluemira.utilities.tools import get_class_from_module + +if TYPE_CHECKING: + from bluemira.base.parameter_frame.typed import ParameterFrameLike + from bluemira.magnets.case_tf import CaseTF + from bluemira.magnets.strand import Strand, SuperconductingStrand + from bluemira.magnets.winding_pack import WindingPack + + +@dataclass +class TFCoilXYDesignerParams(ParameterFrame): + """ + Parameters needed for all aspects of the TF coil design + """ + + # base params + R0: Parameter[float] + """Major radius [m]""" + B0: Parameter[float] + """Magnetic field at R0 [T]""" + A: Parameter[float] + """Aspect ratio""" + n_TF: Parameter[float] + """Number of TF coils""" + ripple: Parameter[float] + """Maximum plasma ripple""" + d: Parameter[float] + """Additional distance to calculate max external radius of inner TF leg""" + S_VV: Parameter[float] + """Vacuum vessel steel limit""" + safety_factor: Parameter[float] + """Allowable stress values""" + B_ref: Parameter[float] + """Reference value for B field (LTS limit) [T]""" + + Iop: Parameter[float] + """Operational current in conductor""" + T_sc: Parameter[float] + """Operational temperature of superconducting cable""" + T_margin: Parameter[float] + """Temperature margin""" + t_delay: Parameter[float] + """Time delay for exponential functions""" + strain: Parameter[float] + """Strain on system""" + + +@dataclass +class StrandParams(ParameterFrame): + """Stand parameters""" + + d_strand: Parameter[float] + operating_temperature: Parameter[float] + + +@dataclass +class CableParams(ParameterFrame): + """Cable parameters""" + + d_cooling_channel: Parameter[float] + void_fraction: Parameter[float] + cos_theta: Parameter[float] + dx: Parameter[float] + + +@dataclass +class CableParamsWithE(CableParams): + """Cable parameters with custom youngs modulus""" + + E: Parameter[float] + + +@dataclass +class SymConductorParams(ParameterFrame): + """Symmetric conductor parameters""" + + dx_ins: Parameter[float] + dx_jacket: Parameter[float] + + +@dataclass +class ConductorParams(SymConductorParams): + """Conductor parameters""" + + dy_ins: Parameter[float] + dy_jacket: Parameter[float] + + +@dataclass +class CaseParams(ParameterFrame): + """Case parameters""" + + # Ri: Parameter[float] + # Rk: Parameter[float] + theta_TF: Parameter[float] + dy_ps: Parameter[float] + dy_vault: Parameter[float] + + +@dataclass +class DerivedTFCoilXYDesignerParams: + """Derived parameters for TF coils cross section design""" + + a: float + """Aspect ratio""" + Ri: float + """External radius of the TF coil case [m].""" + Re: float + B_TF_i: float + pm: float + t_z: float + T_op: float + s_y: float + n_cond: int + min_gap_x: float + I_fun: Callable[[float], float] + B_fun: Callable[[float], float] + + +@dataclass +class OptimisationConfig: + """Optimisation configuration""" + + t0: float + """Initial time""" + Tau_discharge: float + """Characteristic time constant""" + hotspot_target_temperature: float + """Target temperature for hotspot for cable optimisiation""" + layout: str + """Cable layout strategy""" + wp_reduction_factor: float + """Fractional reduction of available toroidal space for Winding Packs""" + n_layers_reduction: int + """Number of layers to remove after each WP""" + bounds_cond_jacket: tuple[float, float] + """Min/max bounds for conductor jacket area optimisation [m²]""" + bounds_dy_vault: tuple[float, float] + """Min/max bounds for the case vault thickness optimisation [m]""" + max_niter: int + """Maximum number of optimisation iterations""" + eps: float + """Convergence threshold for the combined optimisation loop.""" + + +@dataclass +class StabilisingStrandRes: + """Cable opt results""" + + cable: ABCCable + solution: Any + info_text: str + + +@dataclass +class TFCoilXY: + """TFCoil cross section solution""" + + case: CaseTF + cable_soln: StabilisingStrandRes + convergence: npt.NDArray + derived_params: DerivedTFCoilXYDesignerParams + op_config: OptimisationConfig + + def plot_I_B(self, ax, n_steps=300): # noqa: N802 + """Plot current and magnetic field evolution in optimisation""" + time_steps = np.linspace( + self.op_config.t0, self.op_config.Tau_discharge, n_steps + ) + I_values = [self.derived_params.I_fun(t) for t in time_steps] # noqa: N806 + B_values = [self.derived_params.B_fun(t) for t in time_steps] + + ax.plot(time_steps, I_values, "g", label="Current [A]") + ax.set_ylabel("Current [A]", color="g", fontsize=10) + ax.tick_params(axis="y", labelcolor="g", labelsize=9) + ax.grid(visible=True) + + ax_right = ax.twinx() + ax_right.plot(time_steps, B_values, "m--", label="Magnetic field [T]") + ax_right.set_ylabel("Magnetic field [T]", color="m", fontsize=10) + ax_right.tick_params(axis="y", labelcolor="m", labelsize=9) + + # Labels + ax.set_xlabel("Time [s]", fontsize=10) + ax.tick_params(axis="x", labelsize=9) + + # Combined legend for both sides + lines, labels = ax.get_legend_handles_labels() + lines2, labels2 = ax_right.get_legend_handles_labels() + ax.legend(lines + lines2, labels + labels2, loc="best", fontsize=9) + + ax.figure.tight_layout() + + def plot_cable_temperature_evolution(self, ax, n_steps=100): + """Plot temperature evolution in optimisation""" + solution = self.cable_soln.solution + + ax.plot(solution.t, solution.y[0], "r*", label="Simulation points") + time_steps = np.linspace( + self.op_config.t0, self.op_config.Tau_discharge, n_steps + ) + ax.plot(time_steps, solution.sol(time_steps)[0], "b", label="Interpolated curve") + ax.grid(visible=True) + ax.set_ylabel("Temperature [K]", fontsize=10) + ax.set_title("Quench temperature evolution", fontsize=11) + ax.legend(fontsize=9) + + ax.tick_params(axis="y", labelcolor="k", labelsize=9) + + props = {"boxstyle": "round", "facecolor": "white", "alpha": 0.8} + ax.text( + 0.65, + 0.5, + self.cable_soln.info_text, + transform=ax.transAxes, + fontsize=9, + verticalalignment="top", + bbox=props, + ) + ax.figure.tight_layout() + + def plot_summary(self, n_steps, *, show=False): + """Plot summary of optimisation""" # noqa: DOC201 + f, (ax_temp, ax_ib) = plt.subplots(2, 1, figsize=(8, 8), sharex=True) + self.plot_cable_temperature_evolution(ax_temp, n_steps) + self.plot_I_B(ax_ib, n_steps * 3) + if show: + plt.show() + return f + + def plot( + self, + ax: plt.Axes | None = None, + *, + show: bool = False, + homogenised: bool = False, + ) -> plt.Axes: + """Plot the full cross section""" # noqa: DOC201 + return self.case.plot(ax=ax, show=show, homogenised=homogenised) + + def plot_convergence(self, *, show: bool = False): + """ + Plot the evolution of thicknesses and error values over optimisation iterations. + + Raises + ------ + RuntimeError + If no convergence data available + """ + iterations = self.convergence[:, 0] + dy_jacket = self.convergence[:, 1] + dy_vault = self.convergence[:, 2] + err_dy_jacket = self.convergence[:, 3] + err_dy_vault = self.convergence[:, 4] + dy_wp_tot = self.convergence[:, 5] + Ri_minus_Rk = self.convergence[:, 6] # noqa: N806 + + _, axs = plt.subplots(2, 1, figsize=(10, 10), sharex=True) + + # Top subplot: Thicknesses + axs[0].plot(iterations, dy_jacket, marker="o", label="dy_jacket [m]") + axs[0].plot(iterations, dy_vault, marker="s", label="dy_vault [m]") + axs[0].plot(iterations, dy_wp_tot, marker="^", label="dy_wp_tot [m]") + axs[0].plot(iterations, Ri_minus_Rk, marker="v", label="Ri - Rk [m]") + axs[0].set_ylabel("Thickness [m]") + axs[0].set_title("Evolution of Jacket, Vault, and WP Thicknesses") + axs[0].legend() + axs[0].grid(visible=True) + + # Bottom subplot: Errors + axs[1].plot(iterations, err_dy_jacket, marker="o", label="err_dy_jacket") + axs[1].plot(iterations, err_dy_vault, marker="s", label="err_dy_vault") + axs[1].set_ylabel("Relative Error") + axs[1].set_xlabel("Iteration") + axs[1].set_title("Evolution of Errors during Optimisation") + axs[1].set_yscale("log") # Log scale for better visibility if needed + axs[1].legend() + axs[1].grid(visible=True) + + plt.tight_layout() + if show: + plt.show() + + +class TFCoilXYDesigner(Designer[TFCoilXY]): + """ + Handles initialisation of TF Coil XY cross section from the individual parts: + - Strands + - Cable + - Conductor + - Winding Pack + - Casing + + Will output a CaseTF object that allows for the access of all constituent parts + and their properties. + """ + + param_cls: type[TFCoilXYDesignerParams] = TFCoilXYDesignerParams + + def __init__( + self, + params: ParameterFrameLike, + build_config: dict, + ): + super().__init__(params=params, build_config=build_config) + + def _derived_values( + self, op_config: OptimisationConfig + ) -> DerivedTFCoilXYDesignerParams: + # Needed params that are calculated using the base params + R0 = self.params.R0.value # noqa: N806 + n_TF = self.params.n_TF.value + B0 = self.params.B0.value + a = R0 / self.params.A.value + Ri = R0 - a - self.params.d.value # noqa: N806 + Re = (R0 + a) * (1 / self.params.ripple.value) ** (1 / n_TF) # noqa: N806 + B_TF_i = 1.08 * (MU_0_2PI * n_TF * (B0 * R0 / MU_0_2PI / n_TF) / Ri) + t_z = 0.5 * np.log(Re / Ri) * MU_0_4PI * n_TF * (B0 * R0 / MU_0_2PI / n_TF) ** 2 + return DerivedTFCoilXYDesignerParams( + a=a, + Ri=Ri, + Re=Re, + B_TF_i=B_TF_i, + pm=B_TF_i**2 / (2 * MU_0), + t_z=t_z, + T_op=self.params.T_sc.value + self.params.T_margin.value, + s_y=1e9 / self.params.safety_factor.value, + n_cond=int( + self.params.B0.value * R0 / MU_0_2PI / n_TF // self.params.Iop.value + ), + # 2 * thickness of the plate before the WP + min_gap_x=2 * (R0 * 2 / 3 * 1e-2), + I_fun=delayed_exp_func( + self.params.Iop.value, + op_config.Tau_discharge, + self.params.t_delay.value, + ), + B_fun=delayed_exp_func( + B_TF_i, op_config.Tau_discharge, self.params.t_delay.value + ), + ) + + def B_TF_r(self, tf_current: float, r: float): + """ + Compute the magnetic field generated by the TF coils, + including ripple correction. + + Parameters + ---------- + tf_current: + Toroidal field coil current [A]. + r : float + Radial position from the tokamak center [m]. + + Returns + ------- + float + Magnetic field intensity [T]. + """ + return 1.08 * (MU_0_2PI * self.params.n_TF.value * tf_current / r) + + def run(self) -> TFCoilXY: + """ + Run the TF coil XY design problem. + + Returns + ------- + case: + TF case object all parts that make it up. + """ + wp_config = self.build_config.get("winding_pack") + n_wp = int(wp_config.get("sets", 1)) + + optimisation_params = OptimisationConfig( + **self.build_config.get("optimisation_params") + ) + derived_params = self._derived_values(optimisation_params) + + # Only a singular type of cable and conductor currently possible + cable = self.optimise_cable_n_stab_ths( + self._make_cable(wp_i=0), + t0=optimisation_params.t0, + tf=optimisation_params.Tau_discharge, + initial_temperature=derived_params.T_op, + target_temperature=optimisation_params.hotspot_target_temperature, + B_fun=derived_params.B_fun, + I_fun=derived_params.I_fun, + bounds=[1, 10000], + ) + conductor = self._make_conductor(cable.cable, wp_i=0) + winding_pack = [ + self._make_winding_pack(conductor, wp_i, wp_config) for wp_i in range(n_wp) + ] + + # param frame optimisation stuff? + case, convergence_array = self.optimise_jacket_and_vault( + self._make_case(winding_pack, derived_params, optimisation_params), + pm=derived_params.pm, + fz=derived_params.t_z, + op_cond=OperationalConditions( + temperature=derived_params.T_op, + magnetic_field=derived_params.B_TF_i, + strain=self.params.strain.value, + ), + allowable_sigma=derived_params.s_y, + bounds_cond_jacket=optimisation_params.bounds_cond_jacket, + bounds_dy_vault=optimisation_params.bounds_dy_vault, + layout=optimisation_params.layout, + wp_reduction_factor=optimisation_params.wp_reduction_factor, + min_gap_x=derived_params.min_gap_x, + n_layers_reduction=optimisation_params.n_layers_reduction, + max_niter=optimisation_params.max_niter, + eps=optimisation_params.eps, + n_conds=derived_params.n_cond, + ) + return TFCoilXY( + case, cable, convergence_array, derived_params, optimisation_params + ) + + @staticmethod + def optimise_cable_n_stab_ths( + cable: ABCCable, + t0: float, + tf: float, + initial_temperature: float, + target_temperature: float, + B_fun: Callable[[float], float], + I_fun: Callable[[float], float], # noqa: N803 + bounds: np.ndarray | None = None, + ) -> StabilisingStrandRes: + """ + Optimise the number of stabiliser strand in the superconducting cable using a + 0-D hot spot criteria. + + Parameters + ---------- + t0: + Initial time [s]. + tf: + Final time [s]. + initial_temperature: + Temperature [K] at initial time. + target_temperature: + Target temperature [K] at final time. + B_fun : + Magnetic field [T] as a time-dependent function. + I_fun : + Current [A] as a time-dependent function. + bounds: + Lower and upper limits for the number of stabiliser strands. + + Returns + ------- + : + The result of the optimisation process. + + Raises + ------ + ValueError + If the optimisiation process does not converge. + + Notes + ----- + - The number of stabiliser strands in the cable is modified directly. + - Cooling material contribution is neglected when applying the hot spot criteria. + """ + # final_temperature_difference modifies n_stab_strand + # which is used in later calculations of area_stab_region + result = minimize_scalar( + fun=cable.final_temperature_difference, + args=(t0, tf, initial_temperature, target_temperature, B_fun, I_fun), + bounds=bounds, + method=None if bounds is None else "bounded", + ) + + if not result.success: + raise ValueError( + "n_stab optimisation did not converge. Check your input parameters " + "or initial bracket." + ) + + # Here we re-ensure the n_stab_strand to be an integer + cable.n_stab_strand = int(np.ceil(cable.n_stab_strand)) + + solution = cable._temperature_evolution( + t0, tf, initial_temperature, B_fun, I_fun + ) + final_temperature = solution.y[0][-1] + + if final_temperature > target_temperature: + bluemira_error( + f"Final temperature ({final_temperature:.2f} K) exceeds target " + f"temperature " + f"({target_temperature} K) even with maximum n_stab = " + f"{cable.n_stab_strand}." + ) + raise ValueError( + "Optimisation failed to keep final temperature ≤ target. " + "Try increasing the upper bound of n_stab or adjusting cable parameters." + ) + bluemira_print(f"Optimal n_stab: {cable.n_stab_strand}") + bluemira_print( + f"Final temperature with optimal n_stab: {final_temperature:.2f} Kelvin" + ) + + return StabilisingStrandRes( + cable, + solution, + ( + f"Target T: {target_temperature:.2f} K\n" + f"Initial T: {initial_temperature:.2f} K\n" + f"SC Strand: {cable.sc_strand.name}\n" + f"n. sc. strand = {cable.n_sc_strand}\n" + f"Stab. strand = {cable.stab_strand.name}\n" + f"n. stab. strand = {cable.n_stab_strand}\n" + ), + ) + + def optimise_jacket_and_vault( + self, + case: CaseTF, + pm: float, + fz: float, + op_cond: OperationalConditions, + allowable_sigma: float, + bounds_cond_jacket: tuple[float, float] | None = None, + bounds_dy_vault: tuple[float, float] | None = None, + layout: str = "auto", + wp_reduction_factor: float = 0.8, + min_gap_x: float = 0.05, + n_layers_reduction: int = 4, + max_niter: int = 10, + eps: float = 1e-8, + n_conds: int | None = None, + ) -> tuple[CaseTF, npt.NDArray]: + """ + Jointly optimise the conductor jacket and case vault thickness + under electromagnetic loading constraints. + + This method performs an iterative optimisation of: + - The cross-sectional area of the conductor jacket. + - The vault radial thickness of the TF coil casing. + + The optimisation loop continues until the relative change in + jacket area and vault thickness drops below the specified + convergence threshold `eps`, or `max_niter` is reached. + + Parameters + ---------- + pm: + Radial magnetic pressure on the conductor [Pa]. + fz: + Axial electromagnetic force on the winding pack [N]. + op_cond: + Operational conditions including temperature, magnetic field, and strain + at which to calculate the material properties. + allowable_sigma: + Maximum allowable stress for structural material [Pa]. + bounds_cond_jacket: + Min/max bounds for conductor jacket area optimisation [m²]. + bounds_dy_vault: + Min/max bounds for the case vault thickness optimisation [m]. + layout: + Cable layout strategy; "auto" or predefined layout name. + wp_reduction_factor: + Reduction factor applied to WP footprint during conductor rearrangement. + min_gap_x: + Minimum spacing between adjacent conductors [m]. + n_layers_reduction: + Number of conductor layers to remove when reducing WP height. + max_niter: + Maximum number of optimisation iterations. + eps: + Convergence threshold for the combined optimisation loop. + n_conds: + Target total number of conductors in the winding pack. If None, the self + number of conductors is used. + + Returns + ------- + : + Case object + : + convergence array + + Notes + ----- + The function modifies the internal state of `conductor` and `self.dy_vault`. + """ + debug_msg = ["Method optimise_jacket_and_vault"] + + # Initialise convergence array + convergence_array = [] + + if n_conds is None: + n_conds = case.n_conductors + + conductor = case.wps[0].conductor + + case._check_wps(case.wps) + + err_conductor_area_jacket = 10000 * eps + err_dy_vault = 10000 * eps + tot_err = err_dy_vault + err_conductor_area_jacket + + convergence_array.append([ + 0, + conductor.dy_jacket, + case.dy_vault, + err_conductor_area_jacket, + err_dy_vault, + case.dy_wp_tot, + case.geometry.variables.Ri.value - case.geometry.variables.Rk.value, + ]) + + damping_factor = 0.3 + + for i in range(1, max_niter): + if tot_err <= eps: + bluemira_print( + f"Optimisation of jacket and vault reached after " + f"{i - 1} iterations. Total error: {tot_err} < {eps}." + ) + + ax = case.plot(show=False, homogenised=False) + ax.set_title("Case design after optimisation") + plt.show() + break + debug_msg.append(f"Internal optimazion - iteration {i}") + + # Store current values + cond_dx_jacket0 = conductor.dx_jacket + case_dy_vault0 = case.dy_vault + + debug_msg.append( + f"before optimisation: conductor jacket area = {conductor.area_jacket}" + ) + cond_area_jacket0 = conductor.area_jacket + t_z_cable_jacket = ( + fz + * case.area_wps_jacket + / (case.area_case_jacket + case.area_wps_jacket) + / case.n_conductors + ) + self.optimise_jacket_conductor( + conductor, + pm, + t_z_cable_jacket, + op_cond, + allowable_sigma, + bounds_cond_jacket, + ) + debug_msg.extend([ + f"t_z_cable_jacket: {t_z_cable_jacket}", + f"after optimisation: conductor jacket area = {conductor.area_jacket}", + ]) + + conductor.dx_jacket = ( + 1 - damping_factor + ) * cond_dx_jacket0 + damping_factor * conductor.dx_jacket + + err_conductor_area_jacket = ( + abs(conductor.area_jacket - cond_area_jacket0) / cond_area_jacket0 + ) + + case.rearrange_conductors_in_wp( + n_conds, + wp_reduction_factor, + min_gap_x, + n_layers_reduction, + layout=layout, + ) + + debug_msg.append(f"before optimisation: case dy_vault = {case.dy_vault}") + case_dy_vault_result = self.optimise_vault_radial_thickness( + case, + pm=pm, + fz=fz, + op_cond=op_cond, + allowable_sigma=allowable_sigma, + bounds=bounds_dy_vault, + ) + + # case.dy_vault = result.x + # print(f"Optimal dy_vault: {case.dy_vault}") + # print(f"Tresca sigma: {case._tresca_stress(pm, fz, T=T, B=B) / 1e6} MPa") + + case.dy_vault = ( + 1 - damping_factor + ) * case_dy_vault0 + damping_factor * case_dy_vault_result.x + + delta_case_dy_vault = abs(case.dy_vault - case_dy_vault0) + err_dy_vault = delta_case_dy_vault / case.dy_vault + tot_err = err_dy_vault + err_conductor_area_jacket + + debug_msg.append( + f"after optimisation: case dy_vault = {case.dy_vault}\n" + f"err_dy_jacket = {err_conductor_area_jacket}\n " + f"err_dy_vault = {err_dy_vault}\n " + f"tot_err = {tot_err}" + ) + + # Store iteration results in convergence array + convergence_array.append([ + i, + conductor.dy_jacket, + case.dy_vault, + err_conductor_area_jacket, + err_dy_vault, + case.dy_wp_tot, + case.geometry.variables.Ri.value - case.geometry.variables.Rk.value, + ]) + + else: + bluemira_warn( + f"Maximum number of optimisation iterations {max_niter} " + f"reached. A total of {tot_err} > {eps} has been obtained." + ) + + return case, np.array(convergence_array) + + @staticmethod + def optimise_jacket_conductor( + conductor: Conductor, + pressure: float, + f_z: float, + op_cond: OperationalConditions, + allowable_sigma: float, + bounds: tuple[float, float] | None = None, + direction: str = "x", + ): + """ + Optimise the jacket dimension of a conductor based on allowable stress using + the Tresca criterion. + + Parameters + ---------- + pressure: + The pressure applied along the specified direction (Pa). + f_z: + The force applied in the z direction, perpendicular to the conductor + cross-section (N). + op_cond: + Operational conditions including temperature, magnetic field, and strain + at which to calculate the material properties. + allowable_sigma: + The allowable stress (Pa) for the jacket material. + bounds: + Optional bounds for the jacket thickness optimisation (default is None). + direction: + The direction along which the pressure is applied ('x' or 'y'). Default is + 'x'. + + Returns + ------- + : + The result of the optimisation process containing information about the + optimal jacket thickness. + + Raises + ------ + ValueError + If the optimisation process did not converge. + + Notes + ----- + This function uses the Tresca yield criterion to optimise the thickness of the + jacket surrounding the conductor. + This function directly update the conductor's jacket thickness along the x + direction to the optimal value. + """ + debug_msg = ["Method optimise_jacket_conductor:"] + + if direction == "x": + debug_msg.append(f"Previous dx_jacket: {conductor.dx_jacket}") + else: + debug_msg.append(f"Previous dy_jacket: {conductor.dy_jacket}") + + method = "bounded" if bounds is not None else None + + if method == "bounded": + debug_msg.append(f"bounds: {bounds}") + + # sigma difference modifies dx_jacket or dy_jacket which intern is used in + # tresca_sigma_jacket calculation + result = minimize_scalar( + fun=conductor.sigma_difference, + args=(pressure, f_z, op_cond, allowable_sigma), + bounds=bounds, + method=method, + options={"xatol": 1e-4}, + ) + + if not result.success: + raise ValueError("Optimisation of the jacket conductor did not converge.") + if direction == "x": + conductor.dx_jacket = result.x + debug_msg.append(f"Optimal dx_jacket: {conductor.dx_jacket}") + else: + conductor.dy_jacket = result.x + debug_msg.append(f"Optimal dy_jacket: {conductor.dy_jacket}") + debug_msg.append( + f"Averaged sigma in the {direction}-direction: " + f"{conductor._tresca_sigma_jacket(pressure, f_z, op_cond) / 1e6} MPa\n" + f"Allowable stress in the {direction}-direction: {allowable_sigma / 1e6} " + f"MPa." + ) + bluemira_debug("\n".join(debug_msg)) + + return result + + @staticmethod + def optimise_vault_radial_thickness( + case: CaseTF, + pm: float, + fz: float, + op_cond: OperationalConditions, + allowable_sigma: float, + bounds: tuple[float, float] | None = None, + ): + """ + Optimise the vault radial thickness of the case + + Parameters + ---------- + pm: + The magnetic pressure applied along the radial direction (Pa). + f_z: + The force applied in the z direction, perpendicular to the case + cross-section (N). + op_cond: + Operational conditions including temperature, magnetic field, and strain + at which to calculate the material properties. + allowable_sigma: + The allowable stress (Pa) for the jacket material. + bounds: + Optional bounds for the jacket thickness optimisation (default is None). + + Returns + ------- + : + The result of the optimisation process containing information about the + optimal vault thickness. + + Raises + ------ + ValueError + If the optimisation process did not converge. + """ + method = None + if bounds is not None: + method = "bounded" + + # sigma difference modifies dy_vault which is then used in tresca_stress + result = minimize_scalar( + fun=case._sigma_difference, + args=(pm, fz, op_cond, allowable_sigma), + bounds=bounds, + method=method, + options={"xatol": 1e-4}, + ) + + if not result.success: + raise ValueError("dy_vault optimisation did not converge.") + + return result + + def _make_strand( + self, wp_i: int, config: dict[str, str | float], params: StrandParams + ) -> Strand: + cls_name = config["class"] + strand_cls = get_class_from_module( + cls_name, default_module="bluemira.magnets.strand" + ) + material_mix = [] + for m in config.get("materials"): + material_data = m["material"] + if isinstance(material_data, str): + raise TypeError( + "Material data must be a Material instance, not a string - " + "TEMPORARY." + ) + material_mix.append( + MaterialFraction(material=material_data, fraction=m["fraction"]) + ) + + return strand_cls( + materials=material_mix, + d_strand=self._check_iterable(wp_i, params.d_strand.value), + operating_temperature=self._check_iterable( + wp_i, params.operating_temperature.value + ), + name=config.get("name", cls_name.rsplit("::", 1)[-1]), + ) + + @staticmethod + def _check_iterable(wp_i: int, param_val: float | Sequence[float]) -> float: + return param_val[wp_i] if isinstance(param_val, Sequence) else param_val + + def _make_cable_cls( + self, + stab_strand: Strand, + sc_strand: SuperconductingStrand, + wp_i: int, + config: dict[str, float | str], + params: CableParams, + ) -> ABCCable: + cls_name = config["class"] + cable_cls = get_class_from_module( + cls_name, default_module="bluemira.magnets.cable" + ) + extras = {} + if issubclass(cable_cls, RectangularCable): + extras["dx"] = self._check_iterable(wp_i, params.dx.value) + if hasattr(params, "E"): + extras["E"] = self._check_iterable(wp_i, params.E.value) + return cable_cls( + sc_strand=sc_strand, + stab_strand=stab_strand, + n_sc_strand=self._check_iterable(wp_i, config["n_sc_strand"]), + n_stab_strand=self._check_iterable(wp_i, config["n_stab_strand"]), + d_cooling_channel=self._check_iterable(wp_i, params.d_cooling_channel.value), + void_fraction=self._check_iterable(wp_i, params.void_fraction.value), + cos_theta=self._check_iterable(wp_i, params.cos_theta.value), + name=config.get("name", cls_name.rsplit("::", 1)[-1]), + **extras, + ) + + def _make_cable(self, wp_i: int) -> ABCCable: + stab_strand_config = self.build_config.get("stabilising_strand") + sc_strand_config = self.build_config.get("superconducting_strand") + cable_config = self.build_config.get("cable") + + stab_strand_params = StrandParams.from_dict(stab_strand_config.pop("params")) + sc_strand_params = StrandParams.from_dict(sc_strand_config.pop("params")) + cable_params = cable_config.pop("params") + + cable_params = ( + CableParamsWithE.from_dict(cable_params) + if "E" in cable_params + else CableParams.from_dict(cable_params) + ) + + stab_strand = self._make_strand(wp_i, stab_strand_config, stab_strand_params) + sc_strand = self._make_strand(wp_i, sc_strand_config, sc_strand_params) + return self._make_cable_cls( + stab_strand, sc_strand, wp_i, cable_config, cable_params + ) + + def _make_conductor(self, cable: ABCCable, wp_i: int = 0) -> Conductor: + # TODO: current functionality requires conductors are the same for both wps + # in future allow for different conductor objects so can vary cable and + # strands between the sets of the winding pack? + config = self.build_config.get("conductor") + + cls_name = config["class"] + conductor_cls = get_class_from_module( + cls_name, default_module="bluemira.magnets.conductor" + ) + if issubclass(conductor_cls, SymmetricConductor): + params = SymConductorParams.from_dict(config.pop("params")) + else: + params = ConductorParams.from_dict(config.pop("params")) + + return conductor_cls( + cable=cable, + mat_jacket=config["jacket_material"], + mat_ins=config["ins_material"], + dx_jacket=self._check_iterable(wp_i, params.dx_jacket.value), + dx_ins=self._check_iterable(wp_i, params.dx_ins.value), + name=config.get("name", cls_name.rsplit("::", 1)[-1]), + **( + {} + if issubclass(conductor_cls, SymmetricConductor) + else { + "dy_jacket": self._check_iterable(wp_i, params.dy_jacket.value), + "dy_ins": self._check_iterable(wp_i, params.dy_ins.value), + } + ), + ) + + def _make_winding_pack( + self, conductor: Conductor, wp_i: int, config: dict[str, float | str] + ) -> WindingPack: + cls_name = config["class"] + winding_pack_cls = get_class_from_module( + cls_name, default_module="bluemira.magnets.winding_pack" + ) + return winding_pack_cls( + conductor=conductor, + nx=int(self._check_iterable(wp_i, config["nx"])), + ny=int(self._check_iterable(wp_i, config["ny"])), + name="winding_pack", + ) + + def _make_case( + self, + wps: list[WindingPack], + derived_params: DerivedTFCoilXYDesignerParams, + optimisation_params: OptimisationConfig, + ) -> CaseTF: + config = self.build_config.get("case") + params = CaseParams.from_dict(config.get("params")) + + cls_name = config["class"] + case_cls = get_class_from_module( + cls_name, default_module="bluemira.magnets.case_tf" + ) + + case = case_cls( + Ri=derived_params.Ri, + theta_TF=params.theta_TF.value, + dy_ps=params.dy_ps.value, + dy_vault=params.dy_vault.value, + mat_case=config["material"], + wps=wps, + name=config.get("name", cls_name.rsplit("::", 1)[-1]), + ) + + # param frame optimisation stuff? + case.rearrange_conductors_in_wp( + n_conductors=derived_params.n_cond, + wp_reduction_factor=optimisation_params.wp_reduction_factor, + min_gap_x=derived_params.min_gap_x, + n_layers_reduction=optimisation_params.n_layers_reduction, + layout=optimisation_params.layout, + ) + return case diff --git a/bluemira/magnets/utils.py b/bluemira/magnets/utils.py index 597dd5cb3c..b382d4c48b 100644 --- a/bluemira/magnets/utils.py +++ b/bluemira/magnets/utils.py @@ -6,12 +6,14 @@ """Utils for magnets""" +from collections.abc import Callable, Sequence + import numpy as np -def serie_r(arr: list | np.ndarray): +def summation(arr: Sequence[float]) -> float: """ - Compute the serie (as for resistance) + Compute the simple summation of the series Parameters ---------- @@ -21,69 +23,41 @@ def serie_r(arr: list | np.ndarray): Returns ------- - Result: float - """ - return np.sum(arr) - + : + the resulting summation -def parall_r(arr: list | np.ndarray): + Notes + ----- + Y = sum(x1...xn) """ - Compute the parallel (as for resistance) - - Parameters - ---------- - arr: - list or numpy array containing the elements on which the parallel - shall be calculated - - Returns - ------- - Result: float - """ - out = 0 - for i in range(len(arr)): - out += 1 / arr[i] - return out**-1 + return np.sum(arr) -def serie_k(arr: list | np.ndarray): +def reciprocal_summation(arr: Sequence[float]) -> float: """ - Compute the serie (as for spring) + Compute the inverse of the summation of a reciprocal series Parameters ---------- arr: - list or numpy array containing the elements on which the serie - shall be calculated + list or numpy array containing the elements on which the serie shall + be calculated Returns ------- - Result: float - """ - out = 0 - for i in range(len(arr)): - out += 1 / arr[i] - return out**-1 - + : + resulting summation -def parall_k(arr: list | np.ndarray): + Notes + ----- + Y = [sum(1/x1 + 1/x2 + 1/x3 ...)]^-1 """ - Compute the parallel (as for spring) - - Parameters - ---------- - arr: - list or numpy array containing the elements on which the parallel - shall be calculated - - Returns - ------- - Result: float - """ - return np.sum(arr) + return (np.sum(1 / np.array(arr))) ** -1 -def delayed_exp_func(x0: float, tau: float, t_delay: float = 0): +def delayed_exp_func( + x0: float, tau: float, t_delay: float = 0 +) -> Callable[[float], float]: """ Delayed Exponential function @@ -91,16 +65,17 @@ def delayed_exp_func(x0: float, tau: float, t_delay: float = 0): Parameters ---------- - x0: float + x0: initial value - tau: float + tau: characteristic time constant - t_delay: float + t_delay: delay time Returns ------- - A Callable - exponential function + : + An exponential function """ diff --git a/bluemira/magnets/winding_pack.py b/bluemira/magnets/winding_pack.py index 441075787e..abe76704b8 100644 --- a/bluemira/magnets/winding_pack.py +++ b/bluemira/magnets/winding_pack.py @@ -6,66 +6,59 @@ """Winding pack module""" -from typing import Any, ClassVar +from typing import Any import matplotlib.pyplot as plt import numpy as np from matproplib import OperationalConditions -from bluemira.magnets.conductor import Conductor, create_conductor_from_dict -from bluemira.magnets.registry import RegistrableMeta +from bluemira.magnets.conductor import Conductor -# Global registries -WINDINGPACK_REGISTRY = {} - -class WindingPack(metaclass=RegistrableMeta): +class WindingPack: """ Represents a winding pack composed of a grid of conductors. Attributes ---------- - conductor : Conductor + conductor: The base conductor type used in the winding pack. - nx : int + nx: Number of conductors along the x-axis. - ny : int + ny: Number of conductors along the y-axis. """ - _registry_: ClassVar[dict] = WINDINGPACK_REGISTRY - _name_in_registry_: ClassVar[str] = "WindingPack" - def __init__( self, conductor: Conductor, nx: int, ny: int, name: str = "WindingPack" ): """ - Initialize a WindingPack instance. + Initialise a WindingPack instance. Parameters ---------- - conductor : Conductor + conductor: The conductor instance. - nx : int - Number of conductors along the x-direction. - ny : int - Number of conductors along the y-direction. - name : str, optional + nx: + Number of conductors along the x-axis. + ny: + Number of conductors along the y-axis. + name: Name of the winding pack instance. """ self.conductor = conductor - self.nx = int(nx) - self.ny = int(ny) + self.nx = nx + self.ny = ny self.name = name @property def dx(self) -> float: - """Return the total width of the winding pack [m].""" + """Return the width of the winding pack [m].""" return self.conductor.dx * self.nx @property def dy(self) -> float: - """Return the total height of the winding pack [m].""" + """Return the height of the winding pack [m].""" return self.conductor.dy * self.ny @property @@ -89,13 +82,13 @@ def Kx(self, op_cond: OperationalConditions) -> float: # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Stiffness along the x-axis [N/m]. """ return self.conductor.Kx(op_cond) * self.ny / self.nx @@ -106,13 +99,13 @@ def Ky(self, op_cond: OperationalConditions) -> float: # noqa: N802 Parameters ---------- - op_cond: OperationalConditions + op_cond: Operational conditions including temperature, magnetic field, and strain at which to calculate the material property. Returns ------- - float + : Stiffness along the y-axis [N/m]. """ return self.conductor.Ky(op_cond) * self.nx / self.ny @@ -123,28 +116,28 @@ def plot( yc: float = 0, *, show: bool = False, - ax=None, - homogenized: bool = True, - ): + ax: plt.Axes | None = None, + homogenised: bool = True, + ) -> plt.Axes: """ Plot the winding pack geometry. Parameters ---------- - xc : float + xc: Center x-coordinate [m]. - yc : float + yc: Center y-coordinate [m]. - show : bool, optional + show: If True, immediately show the plot. - ax : matplotlib.axes.Axes, optional + ax: Axes object to draw on. - homogenized : bool, optional + homogenised: If True, plot as a single block. Otherwise, plot individual conductors. Returns ------- - matplotlib.axes.Axes + : Axes object containing the plot. """ if ax is None: @@ -164,7 +157,7 @@ def plot( ax.fill(points_ext[:, 0], points_ext[:, 1], "gold", snap=False) ax.plot(points_ext[:, 0], points_ext[:, 1], "k") - if not homogenized: + if not homogenised: for i in range(self.nx): for j in range(self.ny): xc_c = xc - self.dx / 2 + (i + 0.5) * self.conductor.dx @@ -175,116 +168,18 @@ def plot( plt.show() return ax - def to_dict(self) -> dict: + def to_dict(self) -> dict[str, Any]: """ - Serialize the WindingPack to a dictionary. + Serialise the WindingPack to a dictionary. Returns ------- - dict - Serialized dictionary of winding pack attributes. + : + Serialised dictionary of winding pack attributes. """ return { - "name_in_registry": getattr( - self, "_name_in_registry_", self.__class__.__name__ - ), "name": self.name, "conductor": self.conductor.to_dict(), "nx": self.nx, "ny": self.ny, } - - @classmethod - def from_dict( - cls, - windingpack_dict: dict[str, Any], - name: str | None = None, - ) -> "WindingPack": - """ - Deserialize a WindingPack from a dictionary. - - Parameters - ---------- - windingpack_dict : dict - Serialized winding pack dictionary. - name : str - Name for the new instance. If None, attempts to use the 'name' field from - the dictionary. - - Returns - ------- - WindingPack - Reconstructed WindingPack instance. - - Raises - ------ - ValueError - If 'name_in_registry' does not match the expected class. - """ - # Validate name_in_registry - name_in_registry = windingpack_dict.get("name_in_registry") - expected_name_in_registry = getattr(cls, "_name_in_registry_", cls.__name__) - - if name_in_registry != expected_name_in_registry: - raise ValueError( - f"Cannot create {cls.__name__} from dictionary with name_in_registry " - f"'{name_in_registry}'. Expected '{expected_name_in_registry}'." - ) - - # Deserialize conductor - conductor = create_conductor_from_dict( - conductor_dict=windingpack_dict["conductor"], - name=None, - ) - - return cls( - conductor=conductor, - nx=windingpack_dict["nx"], - ny=windingpack_dict["ny"], - name=name or windingpack_dict.get("name"), - ) - - -def create_wp_from_dict( - windingpack_dict: dict[str, Any], - name: str | None = None, -) -> WindingPack: - """ - Factory function to create a WindingPack or its subclass from a serialized - dictionary. - - Parameters - ---------- - windingpack_dict : dict - Dictionary containing serialized winding pack data. - Must include a 'name_in_registry' field matching a registered class. - name : str, optional - Optional name override for the reconstructed WindingPack. - - Returns - ------- - WindingPack - An instance of the appropriate WindingPack subclass. - - Raises - ------ - ValueError - If 'name_in_registry' is missing from the dictionary. - If no matching registered class is found. - """ - name_in_registry = windingpack_dict.get("name_in_registry") - if name_in_registry is None: - raise ValueError( - "Serialized winding pack dictionary must contain a 'name_in_registry' field." - ) - - cls = WINDINGPACK_REGISTRY.get(name_in_registry) - if cls is None: - available = ", ".join(WINDINGPACK_REGISTRY.keys()) - raise ValueError( - f"No registered winding pack class with registration name '" - f"{name_in_registry}'. " - f"Available: {available}." - ) - - return cls.from_dict(windingpack_dict=windingpack_dict, name=name) diff --git a/examples/magnets/example_tf_creation.py b/examples/magnets/example_tf_creation.py new file mode 100644 index 0000000000..f50b1fc5b9 --- /dev/null +++ b/examples/magnets/example_tf_creation.py @@ -0,0 +1,119 @@ +# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza +# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh +# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short +# +# SPDX-License-Identifier: LGPL-2.1-or-later + +""" +Example script demonstrating the design of the TF coil xy cross section. +This involves the design and optimisation of each module: strand, cable, +conductor, winding pack and casing. +""" + +from eurofusion_materials.library.magnet_branch_mats import ( + COPPER_100, + COPPER_300, + DUMMY_INSULATOR_MAG, + NB3SN_MAG, + SS316_LN_MAG, +) + +from bluemira.magnets.tfcoil_designer import TFCoilXYDesigner + +config = { + "stabilising_strand": { + "class": "Strand", + "materials": [{"material": COPPER_300, "fraction": 1.0}], + "params": { + "d_strand": {"value": 1.0e-3, "unit": "m"}, + "operating_temperature": {"value": 5.7, "unit": "K"}, + }, + }, + "superconducting_strand": { + "class": "SuperconductingStrand", + "materials": [ + {"material": NB3SN_MAG, "fraction": 0.5}, + {"material": COPPER_100, "fraction": 0.5}, + ], + "params": { + "d_strand": {"value": 1.0e-3, "unit": "m"}, + "operating_temperature": {"value": 5.7, "unit": "K"}, + }, + }, + "cable": { + "class": "RectangularCable", + "n_sc_strand": 321, + "n_stab_strand": 476, + "params": { + "d_cooling_channel": {"value": 0.01, "unit": "m"}, + "void_fraction": {"value": 0.7, "unit": ""}, + "cos_theta": {"value": 0.97, "unit": ""}, + "dx": {"value": 0.017324217577247843 * 2, "unit": "m"}, + "E": {"value": 0.1e9, "unit": ""}, + }, + }, + "conductor": { + "class": "SymmetricConductor", + "jacket_material": SS316_LN_MAG, + "ins_material": DUMMY_INSULATOR_MAG, + "params": { + "dx_jacket": {"value": 0.0015404278406243683 * 2, "unit": "m"}, + # "dy_jacket": 0.0, + "dx_ins": {"value": 0.0005 * 2, "unit": "m"}, + # "dy_ins": 0.0, + }, + }, + "winding_pack": { + "class": "WindingPack", + "sets": 2, + "nx": [25, 18], + "ny": [6, 1], + }, + "case": { + "class": "TrapezoidalCaseTF", + "material": SS316_LN_MAG, + "params": { + # "Ri": {"value": 3.708571428571428, "unit": "m"}, + # "Rk": {"value": 0, "unit": "m"}, + "theta_TF": {"value": 22.5, "unit": "deg"}, + "dy_ps": {"value": 0.028666666666666667 * 2, "unit": "m"}, + "dy_vault": {"value": 0.22647895819808084 * 2, "unit": "m"}, + }, + }, + "optimisation_params": { + "t0": 0, + "Tau_discharge": 20, + "hotspot_target_temperature": 250, + "layout": "auto", + "wp_reduction_factor": 0.75, + "n_layers_reduction": 4, + "bounds_cond_jacket": (1e-5, 0.2), + "bounds_dy_vault": (0.1, 2), + "max_niter": 100, + "eps": 1e-6, + }, +} + +params = { + # base + "R0": {"value": 8.6, "unit": "m"}, + "B0": {"value": 4.39, "unit": "T"}, + "A": {"value": 2.8, "unit": "dimensionless"}, + "n_TF": {"value": 16, "unit": "dimensionless"}, + "ripple": {"value": 6e-3, "unit": "dimensionless"}, + "d": {"value": 1.82, "unit": "m"}, + "S_VV": {"value": 100e6, "unit": "dimensionless"}, + "safety_factor": {"value": 1.5 * 1.3, "unit": "dimensionless"}, + "B_ref": {"value": 15, "unit": "T"}, + # misc params + "Iop": {"value": 70.0e3, "unit": "A"}, + "T_sc": {"value": 4.2, "unit": "K"}, + "T_margin": {"value": 1.5, "unit": "K"}, + "t_delay": {"value": 3, "unit": "s"}, + "strain": {"value": 0.0055, "unit": ""}, +} + +tf_coil_xy = TFCoilXYDesigner(params=params, build_config=config).execute() +tf_coil_xy.plot(show=True, homogenised=False) +tf_coil_xy.plot_convergence(show=True) +tf_coil_xy.plot_summary(100, show=True) diff --git a/examples/magnets/example_tf_wp_from_dict.py b/examples/magnets/example_tf_wp_from_dict.py deleted file mode 100644 index 2ad3ae77e6..0000000000 --- a/examples/magnets/example_tf_wp_from_dict.py +++ /dev/null @@ -1,261 +0,0 @@ -import matplotlib.pyplot as plt -from bluemira.base.look_and_feel import bluemira_print -import numpy as np -from eurofusion_materials.library.magnet_branch_mats import ( - COPPER_100, - COPPER_300, - DUMMY_INSULATOR_MAG, - NB3SN_MAG, - SS316_LN_MAG, -) -from matproplib import OperationalConditions - -op_cond = OperationalConditions(temperature=5.7, magnetic_field=10.0, strain=0.0055) - -from bluemira.base.constants import MU_0, MU_0_2PI, MU_0_4PI -from bluemira.magnets.case_tf import create_case_tf_from_dict - -case_tf_dict = { - "name_in_registry": "TrapezoidalCaseTF", - "name": "TrapezoidalCaseTF", - "Ri": 3.708571428571428, - "dy_ps": 0.05733333333333333, - "dy_vault": 0.4529579163961617, - "theta_TF": 22.5, - "mat_case": SS316_LN_MAG, - "WPs": [ - { - "name_in_registry": "WindingPack", - "name": "WindingPack", - "conductor": { - "name_in_registry": "SymmetricConductor", - "name": "SymmetricConductor", - "cable": { - "name_in_registry": "DummyRectangularCableLTS", - "name": "DummyRectangularCableLTS", - "n_sc_strand": 321, - "n_stab_strand": 476, - "d_cooling_channel": 0.01, - "void_fraction": 0.7, - "cos_theta": 0.97, - "sc_strand": { - "name_in_registry": "SuperconductingStrand", - "name": "Nb3Sn_strand", - "d_strand": 0.001, - "temperature": 5.7, - "materials": [ - {"material": NB3SN_MAG, "fraction": 0.5}, - {"material": COPPER_100, "fraction": 0.5}, - ], - }, - "stab_strand": { - "name_in_registry": "Strand", - "name": "Stabilizer", - "d_strand": 0.001, - "temperature": 5.7, - "materials": [{"material": COPPER_300, "fraction": 1.0}], - }, - "dx": 0.034648435154495685, - "aspect_ratio": 1.2, - }, - "mat_jacket": SS316_LN_MAG, - "mat_ins": DUMMY_INSULATOR_MAG, - "dx_jacket": 0.0030808556812487366, - "dx_ins": 0.001, - }, - "nx": 25, - "ny": 6, - }, - { - "name_in_registry": "WindingPack", - "name": "WindingPack", - "conductor": { - "name_in_registry": "SymmetricConductor", - "name": "SymmetricConductor", - "cable": { - "name_in_registry": "DummyRectangularCableLTS", - "name": "DummyRectangularCableLTS", - "n_sc_strand": 321, - "n_stab_strand": 476, - "d_cooling_channel": 0.01, - "void_fraction": 0.7, - "cos_theta": 0.97, - "sc_strand": { - "name_in_registry": "SuperconductingStrand", - "name": "Nb3Sn_strand", - "d_strand": 0.001, - "temperature": 5.7, - "materials": [ - {"material": NB3SN_MAG, "fraction": 0.5}, - {"material": COPPER_100, "fraction": 0.5}, - ], - }, - "stab_strand": { - "name_in_registry": "Strand", - "name": "Stabilizer", - "d_strand": 0.001, - "temperature": 5.7, - "materials": [{"material": COPPER_300, "fraction": 1.0}], - }, - "dx": 0.034648435154495685, - "aspect_ratio": 1.2, - }, - "mat_jacket": SS316_LN_MAG, - "mat_ins": DUMMY_INSULATOR_MAG, - "dx_jacket": 0.0030808556812487366, - "dx_ins": 0.001, - }, - "nx": 18, - "ny": 1, - }, - ], -} - -case_tf = create_case_tf_from_dict(case_tf_dict) - -case_tf.plot(show=True, homogenized=False) - -# Machine parameters (should match the original setup) -R0 = 8.6 -B0 = 4.39 -A = 2.8 -n_TF = 16 -ripple = 6e-3 -# operational current per conductor -Iop = 70.0e3 -# Safety factor to be considered on the allowable stress -safety_factor = 1.5 * 1.3 - -# Derived values -a = R0 / A -d = 1.82 -Ri = R0 - a - d -Re = (R0 + a) * (1 / ripple) ** (1 / n_TF) -B_TF_i = 1.08 * (MU_0_2PI * n_TF * (B0 * R0 / MU_0_2PI / n_TF) / Ri) -pm = B_TF_i**2 / (2 * MU_0) -t_z = 0.5 * np.log(Re / Ri) * MU_0_4PI * n_TF * (B0 * R0 / MU_0_2PI / n_TF) ** 2 -T_sc = 4.2 -T_margin = 1.5 -T_op = T_sc + T_margin -S_Y = 1e9 / safety_factor -n_cond = int(np.floor((B0 * R0 / MU_0_2PI / n_TF) / Iop)) - -# Layout and WP parameters -layout = "auto" -wp_reduction_factor = 0.75 -min_gap_x = 2 * (R0 * 2 / 3 * 1e-2) # 2 * dr_plasma_side -n_layers_reduction = 4 - -# Optimization parameters already defined earlier -bounds_cond_jacket = np.array([1e-5, 0.2]) -bounds_dy_vault = np.array([0.1, 2]) -max_niter = 100 -err = 1e-6 - -# optimize number of stabilizer strands -sc_strand = case_tf.WPs[0].conductor.cable.sc_strand -op_cond = OperationalConditions(temperature=T_op, magnetic_field=B_TF_i, strain=0.0055) -Ic_sc = sc_strand.Ic(op_cond) -case_tf.WPs[0].conductor.cable.n_sc_strand = int(np.ceil(Iop / Ic_sc)) - -from bluemira.magnets.utils import delayed_exp_func - -Tau_discharge = 20 # [s] -t_delay = 3 # [s] -t0 = 0 # [s] -hotspot_target_temperature = 250.0 # [K] - -tf = Tau_discharge -T_for_hts = T_op -I_fun = delayed_exp_func(Iop, Tau_discharge, t_delay) -B_fun = delayed_exp_func(B_TF_i, Tau_discharge, t_delay) - - -import time - -t = time.time() -case_tf.WPs[0].conductor.cable.optimize_n_stab_ths( - t0, - tf, - T_for_hts, - hotspot_target_temperature, - B_fun, - I_fun, - bounds=[1, 10000], - show=False, -) - -print("Time taken for optimization:", time.time() - t) - -# Optimize case with structural constraints -case_tf.optimize_jacket_and_vault( - pm=pm, - fz=t_z, - op_cond=OperationalConditions( - temperature=T_op, magnetic_field=B_TF_i, strain=0.0055 - ), - allowable_sigma=S_Y, - bounds_cond_jacket=bounds_cond_jacket, - bounds_dy_vault=bounds_dy_vault, - layout=layout, - wp_reduction_factor=wp_reduction_factor, - min_gap_x=min_gap_x, - n_layers_reduction=n_layers_reduction, - max_niter=max_niter, - eps=err, - n_conds=n_cond, -) - -case_tf.plot_convergence() - -show = True -homogenized = True -if show: - scalex = np.array([2, 1]) - scaley = np.array([1, 1.2]) - - ax = case_tf.plot(homogenized=homogenized) - ax.set_aspect("equal") - - # Fix the x and y limits - ax.set_xlim(-scalex[0] * case_tf.dx_i, scalex[1] * case_tf.dx_i) - ax.set_ylim(scaley[0] * 0, scaley[1] * case_tf.Ri) - - deltax = [-case_tf.dx_i / 2, case_tf.dx_i / 2] - - ax.plot( - [-scalex[0] * case_tf.dx_i, -case_tf.dx_i / 2], [case_tf.Ri, case_tf.Ri], "k:" - ) - - for i in range(len(case_tf.WPs)): - ax.plot( - [-scalex[0] * case_tf.dx_i, -case_tf.dx_i / 2], - [case_tf.R_wp_i[i], case_tf.R_wp_i[i]], - "k:", - ) - - ax.plot( - [-scalex[0] * case_tf.dx_i, -case_tf.dx_i / 2], - [case_tf.R_wp_k[-1], case_tf.R_wp_k[-1]], - "k:", - ) - ax.plot( - [-scalex[0] * case_tf.dx_i, -case_tf.dx_i / 2], [case_tf.Rk, case_tf.Rk], "k:" - ) - - ax.set_title("Equatorial cross section of the TF WP") - ax.set_xlabel("Toroidal direction [m]") - ax.set_ylabel("Radial direction [m]") - - plt.show() - -bluemira_print("Convergence should be: 9.020308301268381e-07 after 11 iterations") - -op_cond = OperationalConditions(temperature=T_op, magnetic_field=B_TF_i, strain=0.0055) -I_sc = case_tf.WPs[0].conductor.cable.sc_strand.Ic(op_cond) -I_max = I_sc * case_tf.WPs[0].conductor.cable.n_sc_strand -I_TF_max = I_max * case_tf.n_conductors -print(I_max) -print(I_TF_max) -I_TF = R0 * B0 / (MU_0_2PI * n_TF) -print(I_TF) diff --git a/examples/magnets/example_tf_wp_optimization.py b/examples/magnets/example_tf_wp_optimization.py deleted file mode 100644 index cdbc0c7bf6..0000000000 --- a/examples/magnets/example_tf_wp_optimization.py +++ /dev/null @@ -1,525 +0,0 @@ -# SPDX-FileCopyrightText: 2021-present M. Coleman, J. Cook, F. Franza -# SPDX-FileCopyrightText: 2021-present I.A. Maione, S. McIntosh -# SPDX-FileCopyrightText: 2021-present J. Morris, D. Short -# -# SPDX-License-Identifier: LGPL-2.1-or-later - -""" -Example script demonstrating the use of the Bluemira WP and TF coil modules -for conductor design, thermal and structural optimization, and case layout visualization. -""" - -# %% md -# # This is an example that shows the application of the WP module for the TF coils -# %% md -# ## Some import -# %% - -import matplotlib.pyplot as plt -import numpy as np -from eurofusion_materials.library.magnet_branch_mats import ( - COPPER_100, - COPPER_300, - DUMMY_INSULATOR_MAG, - NB3SN_MAG, - SS316_LN_MAG, -) -from matplotlib import cm - -# get some materials from EUROfusion materials library -from matproplib import OperationalConditions - -from bluemira.base.constants import MU_0, MU_0_2PI, MU_0_4PI -from bluemira.base.look_and_feel import bluemira_print -from bluemira.magnets.cable import ( - DummyRectangularCableHTS, - DummyRectangularCableLTS, -) -from bluemira.magnets.case_tf import TrapezoidalCaseTF -from bluemira.magnets.conductor import SymmetricConductor -from bluemira.magnets.init_magnets_registry import register_all_magnets -from bluemira.magnets.strand import create_strand_from_dict -from bluemira.magnets.utils import ( - delayed_exp_func, -) -from bluemira.magnets.winding_pack import WindingPack - -# %% -# cache all the magnets classes for future use -register_all_magnets() - - -# %% md -# ## Plot options -# %% -# Enable interactive mode -# %matplotlib notebook - -show = True -homogenized = False -# %% md -# -# ## Input (and derived) values -# *Note: these values shall be provided internally by bluemira code (as reactor -# settings or as information coming from the TF coils builder)

-# -# **Machine (generic)** -# %% -R0 = 8.6 # [m] major machine radius -B0 = 4.39 # [T] magnetic field @R0 -A = 2.8 # machine aspect ratio -n_TF = 16 # number of TF coils -ripple = 6e-3 # requirement on the maximum plasma ripple - -a = R0 / A # minor radius -# %% md -# **Inputs for the TF coils** -# %% -d = 1.82 # additional distance to calculate the max external radius of the inner TF leg -Iop = 70.0e3 # operational current in each conductor -dr_plasma_side = R0 * 2 / 3 * 1e-2 # thickness of the plate before the WP -T_sc = 4.2 # operational temperature of superconducting cable -T_margin = 1.5 # temperature margin -T_op = T_sc + T_margin # temperature considered for the superconducting cable -t_delay = 3 # [s] -t0 = 0 # [s] -hotspot_target_temperature = 250.0 # [K] - -Ri = R0 - a - d # [m] max external radius of the internal TF leg -Re = (R0 + a) * (1 / ripple) ** ( - 1 / n_TF -) # [m] max internal radius of the external TF leg -I_TF = B0 * R0 / MU_0_2PI / n_TF # total current in each TF coil - - -# magnetic field generated by the TF coils with a correction factor that takes into -# account the ripple -def B_TF_r(I_TF, n_TF, r): - """ - Compute the magnetic field generated by the TF coils, including ripple correction. - - Parameters - ---------- - I_TF : float - Toroidal field coil current [A]. - n_TF : int - Number of toroidal field coils. - r : float - Radial position from the tokamak center [m]. - - Returns - ------- - float - Magnetic field intensity [T]. - """ - return 1.08 * (MU_0_2PI * n_TF * I_TF / r) - - -# max magnetic field on the inner TF leg -B_TF_i = B_TF_r(I_TF, n_TF, Ri) -# magnetic pressure on the inner TF leg -pm = B_TF_i**2 / (2 * MU_0) - -# vertical tension acting on the equatorial section of inner TF leg -# i.e. half of the whole F_Z -t_z = 0.5 * np.log(Re / Ri) * MU_0_4PI * n_TF * I_TF**2 - -n_cond = np.floor(I_TF / Iop) # minimum number of conductors -bluemira_print(f"Total number of conductor: {n_cond}") -# %% md -# ***Additional data*** -# %% -R_VV = Ri * 1.05 # Vacuum vessel radius -S_VV = 100e6 # Vacuum vessel steel limit - -# allowable stress values -safety_factor = 1.5 * 1.3 -S_Y = 1e9 / safety_factor # [Pa] steel allowable limit - -# %% md -# ## Calculation of the maximum discharge time for the TF coils -# %% -# inductance (here approximated... better estimation in bluemira) -L = MU_0 * R0 * (n_TF * n_cond) ** 2 * (1 - np.sqrt(1 - (R0 - Ri) / R0)) / n_TF * 1.1 -# Magnetic energy -Wm = 1 / 2 * L * n_TF * Iop**2 * 1e-9 -# Maximum tension... (empirical formula from Lorenzo... find a generic equation) -V_MAX = (7 * R0 - 3) / 6 * 1.1e3 -# Discharge characteristic times -Tau_discharge1 = L * Iop / V_MAX -Tau_discharge2 = B0 * I_TF * n_TF * (R0 / A) ** 2 / (R_VV * S_VV) -# Discharge characteristic time to be considered in the following -Tau_discharge = max([Tau_discharge1, Tau_discharge2]) -tf = Tau_discharge -bluemira_print(f"Maximum TF discharge time: {tf}") -# %% md -# ## Current and magnetic field behaviour during discharge -# -# %% -I_fun = delayed_exp_func(Iop, Tau_discharge, t_delay) -B_fun = delayed_exp_func(B_TF_i, Tau_discharge, t_delay) - -# Create a time array from 0 to 3*Tau_discharge -t = np.linspace(0, 3 * Tau_discharge, 500) -I_data = np.array([I_fun(t_i) for t_i in t]) -B_data = np.array([B_fun(t_i) for t_i in t]) - -# Create a figure and axis -fig, ax1 = plt.subplots() - -# Plot I_fun -ax1.plot(t, I_data, "b-o", label="Current (I)", markevery=30) -ax1.set_xlabel("Time [s]") -ax1.set_ylabel("Current (I) [A]", color="b") -ax1.tick_params("y", colors="b") - -# Create a twin y-axis for the magnetic field B -ax2 = ax1.twinx() -ax2.plot(t, B_data, "r:s", label="Magnetic Field (B)", markevery=30) -ax2.set_ylabel("Magnetic Field (B) [T]", color="r") -ax2.tick_params("y", colors="r") - -# Add grid and title -ax1.grid(visible=True) -plt.title("Current (I) and Magnetic Field (B) vs Time") - -# Show the plot -plt.show() -# %% md -# ### Define materials (at the beginning the conductor is defined with a dummy number -# of stabilizer strands) -# %% -# create the strand -# Define strand configurations as dictionaries -sc_strand_dict = { - "name_in_registry": "SuperconductingStrand", - "name": "Nb3Sn_strand", - "d_strand": 1.0e-3, - "temperature": T_op, - "materials": [ - {"material": NB3SN_MAG, "fraction": 0.5}, - {"material": COPPER_100, "fraction": 0.5}, - ], -} - -stab_strand_dict = { - "name_in_registry": "Strand", - "name": "Stabilizer", - "d_strand": 1.0e-3, - "temperature": T_op, - "materials": [{"material": COPPER_300, "fraction": 1.0}], -} - -# Create strand objects using class-based factory methods -sc_strand = create_strand_from_dict(name="Nb3Sn_strand", strand_dict=sc_strand_dict) -stab_strand = create_strand_from_dict(name="Stabilizer", strand_dict=stab_strand_dict) - -# plot the critical current in a range of B between [10,16] -Bt_arr = np.linspace(10, 16, 100) -sc_strand.plot_Ic_B(Bt_arr, temperature=(T_sc + T_margin)) - -# %% md -# #### Plot number of conductor vs Iop -# %% -# Define the range of operating current (Iop) values -Iop_range = ( - np.linspace(30, 100, 100) * 1e3 -) # 5 equally spaced values between 30 and 100 A -op_conds = [ - OperationalConditions(temperature=T_sc + T_margin, magnetic_field=Bti) - for Bti in Bt_arr -] -Ic_sc_arr = np.array([sc_strand.Ic(op) for op in op_conds]) - -# Create a colormap to assign colors to different Iop values -colors = cm.viridis(np.linspace(0, 1, len(Iop_range))) # Use the 'viridis' colormap - -# Create a figure and axis -fig, ax = plt.subplots() - -# Plot the number of superconducting strands as a function of B for different Iop values -for Iop_ref, color in zip(Iop_range, colors, strict=False): - n_sc_strand = Iop_ref / Ic_sc_arr # Calculate number of strands - ax.plot(Bt_arr, n_sc_strand, color=color, label=f"Iop = {Iop} A") - # Add plot title, axis labels, and grid -ax.set_title("Number of strands as function of B and Iop") # Title -ax.set_xlabel("B [T]") # X-axis label -ax.set_ylabel("Nc strands") # Y-axis label -ax.grid(visible=True) - -# Create a ScalarMappable to map colors to the colorbar -sm = plt.cm.ScalarMappable( - cmap="viridis", norm=plt.Normalize(vmin=Iop_range.min(), vmax=Iop_range.max()) -) -sm.set_array([]) # Dummy array for the ScalarMappable - -# Add the colorbar to the figure -cbar = fig.colorbar(sm, ax=ax) -cbar.set_label("Iop [A]") # Label the colorbar - -# Show the plot -plt.show() - -# %% md -# **Calculate number of superconducting strands considering the strand critical -# current at B_TF_i and T_sc + T_margin** -# %% -op_cond = OperationalConditions(temperature=T_op, magnetic_field=B_TF_i) -Ic_sc = sc_strand.Ic(op_cond) -n_sc_strand = int(np.ceil(Iop / Ic_sc)) - -########################################################### -dx = 0.05 # cable length... just a dummy value -B_ref = 15 # [T] Reference B field value (limit for LTS) - -if B_TF_i < B_ref: - cable = DummyRectangularCableLTS( - dx=dx, - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=n_sc_strand, - n_stab_strand=500, - d_cooling_channel=1e-2, - void_fraction=0.7, - cos_theta=0.97, - ) -else: - cable = DummyRectangularCableHTS( - dx=dx, - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=n_sc_strand, - n_stab_strand=500, - d_cooling_channel=1e-2, - void_fraction=0.7, - cos_theta=0.97, - ) -cable.plot(0, 0, show=True) -bluemira_print(f"cable area: {cable.area}") -# %% md -# -# %% md -# ***Change cable aspect ratio*** -# %% -aspect_ratio = 1.2 -cable.set_aspect_ratio( - aspect_ratio -) # This adjusts the cable dimensions while maintaining the total cross-sectional area. -cable.plot(0, 0, show=True) -bluemira_print(f"cable area: {cable.area}") - -# operational_point = {"temperature": 5.7, "B": B(0)} - -mats = [NB3SN_MAG, COPPER_100, COPPER_300, SS316_LN_MAG] -mat_names = ["nb3sn", "copper100", "copper300", "ss316"] -temperatures = np.linspace(5, 250, 500) -magnetic_field = B_fun(0) - -# Prepare plots -# Adjusted for 5 plots (3x2 grid) -fig, axes = plt.subplots(3, 2, figsize=(12, 15)) -fig.suptitle("Material Properties vs Temperature", fontsize=16) - -ax1, ax2, ax3, ax4, ax5 = axes.flatten()[:5] # Extracting the first five axes - -for mat, name in zip(mats, mat_names, strict=False): - # Calculate properties over the temperature range - op_conds = [ - OperationalConditions(temperature=T, magnetic_field=magnetic_field) - for T in temperatures - ] - density = np.array([mat.density(op) for op in op_conds]) - cp_per_mass = np.array([mat.specific_heat_capacity(op) for op in op_conds]) - cp_per_volume = cp_per_mass * density - erho = np.array([mat.electrical_resistivity(op) for op in op_conds]) - E = np.array([mat.youngs_modulus(op) for op in op_conds]) - - # Plot density - ax1.plot(temperatures, density, label=name) - # Plot Cp [J/Kg/K] - ax2.plot(temperatures, cp_per_mass, label=name) - # Plot Cp [J/K/m³] - ax3.plot(temperatures, cp_per_volume, label=name) - # Plot erho [Ohm m] - ax4.plot(temperatures, erho, label=name) - # Plot E (assuming Young's modulus) - ax5.plot(temperatures, E, label=name) - -# Configure plots -ax1.set_title("Density [kg/m³] vs Temperature") -ax1.set_xlabel("Temperature [K]") -ax1.set_ylabel("Density [kg/m³]") -ax1.legend() -ax1.grid(visible=True) - -ax2.set_title("Heat Capacity per Mass [J/Kg/K] vs Temperature") -ax2.set_xlabel("Temperature [K]") -ax2.set_ylabel("Cp [J/Kg/K]") -ax2.legend() -ax2.grid(visible=True) - -ax3.set_title("Heat Capacity per Volume [J/K/m³] vs Temperature") -ax3.set_xlabel("Temperature [K]") -ax3.set_ylabel("Cp [J/K/m³]") -ax3.legend() -ax3.grid(visible=True) - -ax4.set_title("Electrical Resistivity [Ohm·m] vs Temperature") -ax4.set_xlabel("Temperature [K]") -ax4.set_ylabel("Resistivity [Ohm·m]") -ax4.legend() -ax4.grid(visible=True) - -ax5.set_title("E vs Temperature") # Modify title according to the meaning of E -ax5.set_xlabel("Temperature [K]") -ax5.set_ylabel("E [Unit]") # Change to correct unit -ax5.legend() -ax5.grid(visible=True) - -plt.tight_layout(rect=[0, 0, 1, 0.96]) -plt.show() - -# %% -# optimize the number of stabilizer strands using the hot spot criteria. -# Note: optimize_n_stab_ths changes adjust cable.dy while maintaining cable.dx. It -# could be possible to add a parameter maintain constant the aspect ratio if -# necessary. -bluemira_print( - f"before optimization: dx_cable = {cable.dx}, aspect ratio = {cable.aspect_ratio}" -) -T_for_hts = T_op -result = cable.optimize_n_stab_ths( - t0, - tf, - T_for_hts, - hotspot_target_temperature, - B_fun, - I_fun, - bounds=[1, 10000], - show=show, -) -bluemira_print( - f"after optimization: dx_cable = {cable.dx}, aspect ratio = {cable.aspect_ratio}" -) - -cable.set_aspect_ratio(aspect_ratio) -bluemira_print( - f"Adjust aspect ratio: dx_cable = {cable.dx}, aspect ratio = {cable.aspect_ratio}" -) - -# %% -########################################################### -# Create a conductor with the specified cable -conductor = SymmetricConductor( - cable=cable, - mat_jacket=SS316_LN_MAG, - mat_ins=DUMMY_INSULATOR_MAG, - dx_jacket=0.01, - dx_ins=1e-3, -) - -# %% -# case parameters -layout = "auto" # "layer" or "pancake" -wp_reduction_factor = 0.75 -min_gap_x = 2 * dr_plasma_side -n_layers_reduction = 4 - -# creation of the case -wp1 = WindingPack(conductor, 1, 1, name=None) # just a dummy WP to create the case - -case = TrapezoidalCaseTF( - Ri=Ri, - dy_ps=dr_plasma_side, - dy_vault=0.7, - theta_TF=360 / n_TF, - mat_case=SS316_LN_MAG, - WPs=[wp1], -) - -print(f"pre-wp reduction factor: {wp_reduction_factor}") - -# arrangement of conductors into the winding pack and case -case.rearrange_conductors_in_wp( - n_conductors=n_cond, - wp_reduction_factor=wp_reduction_factor, - min_gap_x=min_gap_x, - n_layers_reduction=n_layers_reduction, - layout=layout, -) - -ax = case.plot(show=False, homogenized=False) -ax.set_title("Case design before optimization") -plt.show() - -bluemira_print(f"Previous number of conductors: {n_cond}") -bluemira_print(f"New number of conductors: {case.n_conductors}") - -# %% md -# ## Optimize cable jacket and case vault thickness -# %% -# Optimization parameters -bounds_cond_jacket = np.array([1e-5, 0.2]) -bounds_dy_vault = np.array([0.1, 2]) -max_niter = 100 -err = 1e-6 - -# case.optimize_jacket_and_vault( -case.optimize_jacket_and_vault( - pm=pm, - fz=t_z, - op_cond=OperationalConditions(temperature=T_op, magnetic_field=B_TF_i), - allowable_sigma=S_Y, - bounds_cond_jacket=bounds_cond_jacket, - bounds_dy_vault=bounds_dy_vault, - layout=layout, - wp_reduction_factor=wp_reduction_factor, - min_gap_x=min_gap_x, - n_layers_reduction=n_layers_reduction, - max_niter=max_niter, - eps=err, - n_conds=n_cond, -) - -if show: - scalex = np.array([2, 1]) - scaley = np.array([1, 1.2]) - - ax = case.plot(homogenized=homogenized) - ax.set_aspect("equal") - - # Fix the x and y limits - ax.set_xlim(-scalex[0] * case.dx_i, scalex[1] * case.dx_i) - ax.set_ylim(scaley[0] * 0, scaley[1] * case.Ri) - - deltax = [-case.dx_i / 2, case.dx_i / 2] - - ax.plot([-scalex[0] * case.dx_i, -case.dx_i / 2], [case.Ri, case.Ri], "k:") - - for i in range(len(case.WPs)): - ax.plot( - [-scalex[0] * case.dx_i, -case.dx_i / 2], - [case.R_wp_i[i], case.R_wp_i[i]], - "k:", - ) - - ax.plot( - [-scalex[0] * case.dx_i, -case.dx_i / 2], - [case.R_wp_k[-1], case.R_wp_k[-1]], - "k:", - ) - ax.plot([-scalex[0] * case.dx_i, -case.dx_i / 2], [case.Rk, case.Rk], "k:") - - ax.set_title("Equatorial cross section of the TF WP") - ax.set_xlabel("Toroidal direction [m]") - ax.set_ylabel("Radial direction [m]") - - plt.show() - - -bluemira_print("Convergence should be: 9.066682976310327e-07 after 68 iterations") -# %% -# new operational current -bluemira_print(f"Operational current after optimization: {I_TF / case.n_conductors}") - -case.plot_convergence() diff --git a/tests/magnets/test_cable.py b/tests/magnets/test_cable.py index 0dd05dada7..fae6344923 100644 --- a/tests/magnets/test_cable.py +++ b/tests/magnets/test_cable.py @@ -5,7 +5,6 @@ # SPDX-License-Identifier: LGPL-2.1-or-later -import matplotlib.pyplot as plt import numpy as np import pytest from eurofusion_materials.library.magnet_branch_mats import ( @@ -15,18 +14,13 @@ from matproplib import OperationalConditions from matproplib.material import MaterialFraction -from bluemira.magnets.cable import ( - DummyRoundCableLTS, - DummySquareCableLTS, - RectangularCable, -) +from bluemira.magnets.cable import RectangularCable, RoundCable, SquareCable from bluemira.magnets.strand import Strand, SuperconductingStrand +from bluemira.magnets.tfcoil_designer import TFCoilXYDesigner DummySteel = SS316_LN_MAG DummySuperconductor = NB3SN_MAG -# -- Pytest Fixtures ---------------------------------------------------------- - @pytest.fixture def sc_strand(): @@ -34,6 +28,7 @@ def sc_strand(): name="SC", materials=[MaterialFraction(material=DummySuperconductor, fraction=1.0)], d_strand=0.001, + operating_temperature=5.7, ) @@ -43,6 +38,7 @@ def stab_strand(): name="Stab", materials=[MaterialFraction(material=DummySteel, fraction=1.0)], d_strand=0.001, + operating_temperature=5.7, ) @@ -55,20 +51,19 @@ def cable(sc_strand, stab_strand): n_sc_strand=10, n_stab_strand=20, d_cooling_channel=0.001, + void_fraction=0.725, + cos_theta=0.97, ) -# -- Core Cable Tests --------------------------------------------------------- - - def test_geometry_and_area(cable): assert cable.dx > 0 assert cable.dy > 0 assert cable.area > 0 assert cable.aspect_ratio > 0 - assert cable.area_cc > 0 - assert cable.area_stab > 0 - assert cable.area_sc > 0 + assert cable.area_cooling_channel > 0 + assert cable.area_stab_region > 0 + assert cable.area_sc_region > 0 def test_material_properties(cable): @@ -86,8 +81,7 @@ def test_str_output(cable): assert "stab strand" in summary -def test_plot(monkeypatch, cable): - monkeypatch.setattr(plt, "show", lambda: None) +def test_plot(cable): ax = cable.plot(show=True) assert hasattr(ax, "fill") @@ -103,7 +97,7 @@ def I_fun(t): # noqa: ARG001 assert result.success -def test_optimize_n_stab_ths(monkeypatch, sc_strand, stab_strand): +def test_optimise_n_stab_ths(sc_strand, stab_strand): cable = RectangularCable( dx=0.01, sc_strand=sc_strand, @@ -111,8 +105,9 @@ def test_optimize_n_stab_ths(monkeypatch, sc_strand, stab_strand): n_sc_strand=10, n_stab_strand=5, d_cooling_channel=0.001, + void_fraction=0.725, + cos_theta=0.97, ) - monkeypatch.setattr(plt, "show", lambda: None) def B_fun(t): # noqa: ARG001 return 5 @@ -120,7 +115,8 @@ def B_fun(t): # noqa: ARG001 def I_fun(t): # noqa: ARG001 return 1000 - result = cable.optimize_n_stab_ths( + result = TFCoilXYDesigner.optimise_cable_n_stab_ths( + cable, t0=0, tf=0.1, initial_temperature=20, @@ -129,49 +125,29 @@ def I_fun(t): # noqa: ARG001 I_fun=I_fun, bounds=(1, 100), ) - assert result.success - - -def test_invalid_parameters(sc_strand, stab_strand): - with pytest.raises(ValueError, match="dx must be positive"): - RectangularCable( - dx=-0.01, - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=10, - n_stab_strand=5, - d_cooling_channel=0.001, - ) - - with pytest.raises(ValueError, match="void_fraction must be between 0 and 1"): - RectangularCable( - dx=0.01, - sc_strand=sc_strand, - stab_strand=stab_strand, - n_sc_strand=10, - n_stab_strand=5, - d_cooling_channel=0.001, - void_fraction=1.5, - ) - - -# -- Square & Round Cable Types ----------------------------------------------- + assert result.solution.success def test_square_and_round_cables(sc_strand, stab_strand): - square = DummySquareCableLTS( + square = SquareCable( sc_strand=sc_strand, stab_strand=stab_strand, n_sc_strand=5, n_stab_strand=5, d_cooling_channel=0.001, + void_fraction=0.725, + cos_theta=0.97, + E=0.1e9, ) - round_ = DummyRoundCableLTS( + round_ = RoundCable( sc_strand=sc_strand, stab_strand=stab_strand, n_sc_strand=5, n_stab_strand=5, d_cooling_channel=0.001, + void_fraction=0.725, + cos_theta=0.97, + E=0.1e9, ) dummy_op_cond = OperationalConditions(temperature=4.0) assert square.dx > 0 @@ -191,7 +167,7 @@ def test_square_and_round_cables(sc_strand, stab_strand): assert np.isclose(round_.Ky(dummy_op_cond), round_.E(dummy_op_cond), rtol=1e-8) -def test_cable_to_from_dict(sc_strand, stab_strand): +def test_cable_to_dict(sc_strand, stab_strand): # Create a RectangularCable for testing cable_original = RectangularCable( dx=0.01, @@ -200,24 +176,11 @@ def test_cable_to_from_dict(sc_strand, stab_strand): n_sc_strand=10, n_stab_strand=5, d_cooling_channel=0.001, + void_fraction=0.725, + cos_theta=0.97, ) # Convert to dictionary - cable_dict = cable_original.to_dict() + cable_dict = cable_original.to_dict(OperationalConditions(temperature=5)) - # Reconstruct from dictionary - cable_reconstructed = RectangularCable.from_dict(cable_dict) - - # Verify key attributes match - assert cable_original.n_sc_strand == cable_reconstructed.n_sc_strand - assert cable_original.n_stab_strand == cable_reconstructed.n_stab_strand - assert cable_original.dx == pytest.approx(cable_reconstructed.dx) - assert cable_original.d_cooling_channel == pytest.approx( - cable_reconstructed.d_cooling_channel - ) - assert cable_original.void_fraction == pytest.approx( - cable_reconstructed.void_fraction - ) - assert cable_original.cos_theta == pytest.approx(cable_reconstructed.cos_theta) - assert cable_original.sc_strand.name == cable_reconstructed.sc_strand.name - assert cable_original.stab_strand.name == cable_reconstructed.stab_strand.name + assert cable_dict["n_sc_strand"] == 10 diff --git a/tests/magnets/test_conductor.py b/tests/magnets/test_conductor.py index c3f1d2c3a9..ccaae0d6fa 100644 --- a/tests/magnets/test_conductor.py +++ b/tests/magnets/test_conductor.py @@ -19,10 +19,6 @@ from bluemira.magnets.conductor import Conductor from bluemira.magnets.strand import Strand, SuperconductingStrand -# --------------------------------------------------------------------------- -# Fixtures -# --------------------------------------------------------------------------- - @pytest.fixture def mat_jacket(): @@ -42,6 +38,7 @@ def sc_strand(): name="SC", materials=[MaterialFraction(material=sc, fraction=1.0)], d_strand=0.001, + operating_temperature=5.7, ) @@ -53,6 +50,7 @@ def stab_strand(): name="Stab", materials=[MaterialFraction(material=stab, fraction=1.0)], d_strand=0.001, + operating_temperature=5.7, ) @@ -65,6 +63,8 @@ def rectangular_cable(sc_strand, stab_strand): n_sc_strand=10, n_stab_strand=10, d_cooling_channel=0.001, + void_fraction=0.725, + cos_theta=0.97, ) @@ -82,11 +82,6 @@ def conductor(rectangular_cable, mat_jacket, mat_ins): ) -# --------------------------------------------------------------------------- -# Tests -# --------------------------------------------------------------------------- - - def test_geometry_and_area(conductor): assert conductor.dx > 0 assert conductor.dy > 0 @@ -105,18 +100,3 @@ def test_plot(monkeypatch, conductor): monkeypatch.setattr(plt, "show", lambda: None) ax = conductor.plot(show=True) assert hasattr(ax, "fill") - - -def test_to_from_dict(conductor): - config = conductor.to_dict() - restored = Conductor.from_dict(config) - - assert restored.name == conductor.name - assert restored.dx_jacket == pytest.approx(conductor.dx_jacket) - assert restored.dy_jacket == pytest.approx(conductor.dy_jacket) - assert restored.dx_ins == pytest.approx(conductor.dx_ins) - assert restored.dy_ins == pytest.approx(conductor.dy_ins) - assert restored.mat_jacket.name == conductor.mat_jacket.name - assert restored.mat_ins.name == conductor.mat_ins.name - assert restored.cable.n_sc_strand == conductor.cable.n_sc_strand - assert restored.cable.sc_strand.name == conductor.cable.sc_strand.name diff --git a/tests/magnets/test_strand.py b/tests/magnets/test_strand.py index af5190fb34..864444e81e 100644 --- a/tests/magnets/test_strand.py +++ b/tests/magnets/test_strand.py @@ -26,34 +26,36 @@ def test_strand_area(): mat = MaterialFraction(material=DummySuperconductor1, fraction=1.0) - strand = Strand(name="test_strand", materials=[mat], d_strand=0.001) + strand = Strand( + name="test_strand", materials=[mat], d_strand=0.001, operating_temperature=5.7 + ) expected_area = np.pi * (0.001**2) / 4 assert np.isclose(strand.area, expected_area) -def test_strand_invalid_diameter(): - mat = MaterialFraction(material=DummySuperconductor1, fraction=1.0) - with pytest.raises(ValueError, match="positive"): - Strand(name="invalid_strand", materials=[mat], d_strand=-0.001) - - def test_superconducting_strand_invalid_materials(): # Two superconductors — should raise ValueError mat1 = MaterialFraction(material=DummySuperconductor1, fraction=0.5) mat2 = MaterialFraction(material=DummySuperconductor2, fraction=0.5) with pytest.raises(ValueError, match="Only one superconductor material"): - SuperconductingStrand(name="invalid", materials=[mat1, mat2]) + SuperconductingStrand( + name="invalid", materials=[mat1, mat2], d_strand=0, operating_temperature=0 + ) # No superconductors — should raise ValueError mat3 = MaterialFraction(material=DummySteel, fraction=1.0) with pytest.raises(ValueError, match="No superconducting material"): - SuperconductingStrand(name="invalid", materials=[mat3]) + SuperconductingStrand( + name="invalid", materials=[mat3], d_strand=0, operating_temperature=0 + ) def test_strand_material_properties(): sc = DummySuperconductor1 mat = MaterialFraction(material=sc, fraction=1.0) - strand = Strand(name="mat_test", materials=[mat], d_strand=0.001) + strand = Strand( + name="mat_test", materials=[mat], d_strand=0.001, operating_temperature=5.7 + ) temperature = 20 op_cond = OperationalConditions(temperature=20)