From 2559815dfe32417b4bfdd0432825db78dd7d1e38 Mon Sep 17 00:00:00 2001 From: james Date: Fri, 13 Sep 2024 07:58:55 +0100 Subject: [PATCH 1/5] =?UTF-8?q?=E2=9C=A8=20Coil=20resistance?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/equilibria/coils/_coil.py | 16 ++++++++++++++++ bluemira/equilibria/coils/_grouping.py | 10 ++++++++++ 2 files changed, 26 insertions(+) diff --git a/bluemira/equilibria/coils/_coil.py b/bluemira/equilibria/coils/_coil.py index e028f55fdd..542db29512 100644 --- a/bluemira/equilibria/coils/_coil.py +++ b/bluemira/equilibria/coils/_coil.py @@ -127,6 +127,7 @@ class Coil(CoilFieldsMixin): "_flag_sizefix", "_j_max", "_number", + "_resistance", "_x", "_x_boundary", "_z", @@ -148,6 +149,7 @@ def __init__( b_max: float = np.nan, discretisation: float = np.nan, n_turns: int = 1, + resistance: float = 0, *, psi_analytic: bool = False, Bx_analytic: bool = True, @@ -172,6 +174,7 @@ def __init__( self.ctype = ctype self.name = name self.n_turns = n_turns + self.resistance = resistance self._number = CoilNumber.generate(self.ctype) if self.name is None: @@ -196,6 +199,7 @@ def __repr__(self): f"{type(self).__name__}({self.name} ctype={self.ctype.name} x={self.x:.2g}" f" z={self.z:.2g} dx={self.dx:.2g} dz={self.dz:.2g}" f" current={self.current:.2g} j_max={self.j_max:.2g} b_max={self.b_max:.2g}" + f" resistance={self.resistance:.2g}" f" discretisation={self.discretisation:.2g})" ) @@ -410,6 +414,16 @@ def b_max(self, value: float): """Set coil max field""" self._b_max = floatify(value) + @property + def resistance(self): + """Get coil resistance""" + return self._resistance + + @resistance.setter + def resistance(self, value: float): + """Set coil resistance""" + self._resistance = value + @discretisation.setter def discretisation(self, value: float): """Set coil discretisation""" @@ -420,6 +434,7 @@ def assign_material( self, j_max: float = NBTI_J_MAX, b_max: float = NBTI_B_MAX, + resistance: float = 0, ) -> None: """ Assigns EM material properties to coil @@ -438,6 +453,7 @@ def assign_material( """ self.j_max = j_max self.b_max = b_max + self.resistance = resistance def get_max_current(self): """Get max current""" diff --git a/bluemira/equilibria/coils/_grouping.py b/bluemira/equilibria/coils/_grouping.py index 047bc0dc18..cc44ace2be 100644 --- a/bluemira/equilibria/coils/_grouping.py +++ b/bluemira/equilibria/coils/_grouping.py @@ -698,6 +698,11 @@ def b_max(self) -> np.ndarray: """Get coil max field""" return self.__getter("b_max") + @property + def resistance(self) -> np.ndarray: + """Get coil resistance""" + return self.__getter("resistance") + @property def discretisation(self) -> np.ndarray: """Get coil discretisations""" @@ -829,6 +834,11 @@ def b_max(self, values: float | Iterable[float]): """Set coil max fields""" self.__setter("b_max", values) + @resistance.setter + def resistance(self, values: float | Iterable[float]): + """Set coil resistance""" + self.__setter("b_max", values) + @discretisation.setter def discretisation(self, values: float | Iterable[float]): """Set coil discretisations""" From 8793af8b097cfcfb7a750cd8cd8e8e98d54e65e5 Mon Sep 17 00:00:00 2001 From: james Date: Fri, 13 Sep 2024 08:03:36 +0100 Subject: [PATCH 2/5] =?UTF-8?q?=E2=9C=A8=20Differential=20greens=20functio?= =?UTF-8?q?ns?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/equilibria/coils/_field.py | 112 +++++++++++++++++- bluemira/equilibria/equilibrium.py | 91 ++++++++++++++- bluemira/equilibria/plasma.py | 175 ++++++++++++++++++++++++++-- bluemira/magnetostatics/greens.py | 164 +++++++++++++++++++++++++- 4 files changed, 523 insertions(+), 19 deletions(-) diff --git a/bluemira/equilibria/coils/_field.py b/bluemira/equilibria/coils/_field.py index 283402aacd..19efc252d4 100644 --- a/bluemira/equilibria/coils/_field.py +++ b/bluemira/equilibria/coils/_field.py @@ -16,7 +16,12 @@ from bluemira.base.constants import MU_0 from bluemira.equilibria.constants import X_TOLERANCE -from bluemira.magnetostatics.greens import greens_Bx, greens_Bz, greens_psi +from bluemira.magnetostatics.greens import ( + greens_Bx, + greens_Bz, + greens_dbz_dx, + greens_psi, +) from bluemira.magnetostatics.semianalytic_2d import ( semianalytic_Bx, semianalytic_Bz, @@ -59,6 +64,50 @@ def __init__( self._Bx_analytic = Bx_analytic self._Bz_analytic = Bz_analytic + def dB_d(self, x: float | np.ndarray, z: float | np.ndarray): + """ + dB_d of Coilset + + Parameters + ---------- + x: + The x values at which to calculate the dB_d response + z: + The z values at which to calculate the dB_d response + sum_coils: + sum over coils + control: + operations on control coils only + + Returns + ------- + : + Differential of magnetic field + """ + return self.dB_d_response(x, z) * self.current + + def dB_d_response(self, x: float | np.ndarray, z: float | np.ndarray): + """ + Unit dB_d of Coilset + + Parameters + ---------- + x: + The x values at which to calculate the dB_d response + z: + The z values at which to calculate the dB_d response + sum_coils: + sum over coils + control: + operations on control coils only + + Returns + ------- + : + Differential of magnetic field response + """ + return self._mix_control_method(x, z, greens_dbz_dx, None, disable_analytic=True) + def psi(self, x: float | np.ndarray, z: float | np.ndarray): """ Calculate poloidal flux at (x, z) @@ -593,6 +642,34 @@ def Bz( """ return self._sum(super().Bz(x, z), sum_coils=sum_coils, control=control) + def dB_d( + self, + x: np.ndarray, + z: np.ndarray, + *, + sum_coils: bool = True, + control: bool = False, + ) -> np.ndarray: + """ + dB_d of Coilset + + Parameters + ---------- + x: + The x values at which to calculate the dB_d response + z: + The z values at which to calculate the dB_d response + sum_coils: + sum over coils + control: + operations on control coils only + + Returns + ------- + Differential of magnetic field + """ + return self._sum(super().dB_d(x, z), sum_coils=sum_coils, control=control) + def psi_response( self, x: np.ndarray, @@ -660,7 +737,7 @@ def Bz_response( control: bool = False, ) -> np.ndarray: """ - Bz of Coilset + Unit Bz of Coilset Parameters ---------- @@ -679,6 +756,37 @@ def Bz_response( """ return self._sum(super().Bz_response(x, z), sum_coils=sum_coils, control=control) + def dB_d_response( + self, + x: np.ndarray, + z: np.ndarray, + *, + sum_coils: bool = False, + control: bool = False, + ) -> np.ndarray: + """ + Unit dB_d of Coilset + + Parameters + ---------- + x: + The x values at which to calculate the dB_d response + z: + The z values at which to calculate the dB_d response + sum_coils: + sum over coils + control: + operations on control coils only + + Returns + ------- + : + Differential of magnetic field response + """ + return self._sum( + super().dB_d_response(x, z), sum_coils=sum_coils, control=control + ) + def control_F(self, coil_grp: CoilGroup, *, control: bool = False) -> np.ndarray: """ Returns the Green's matrix element for the coil mutual force. diff --git a/bluemira/equilibria/equilibrium.py b/bluemira/equilibria/equilibrium.py index 967b95053e..55f84fe26d 100644 --- a/bluemira/equilibria/equilibrium.py +++ b/bluemira/equilibria/equilibrium.py @@ -524,6 +524,7 @@ def _remap_greens(self): self._psi_green = self.coilset.psi_response(self.x, self.z) self._bx_green = self.coilset.Bx_response(self.x, self.z) self._bz_green = self.coilset.Bz_response(self.x, self.z) + self._db_green = self.coilset.dB_d_response(self.x, self.z) def get_coil_forces(self) -> npt.NDArray[np.float64]: """ @@ -836,6 +837,54 @@ def Bp( return np.hypot(self.Bx(x, z), self.Bz(x, z)) + def dBx_dz(self, x: npt.ArrayLike, z: npt.ArrayLike): + """ + Total differential of the radial magnetic field at point (x, z) from coils + + Parameters + ---------- + x: + Radial coordinates for which to return Bz. If None, returns values + at all grid points + z: + Vertical coordinates for which to return Bz. If None, returns values + at all grid points + + Returns + ------- + : + Differential of the radial magnetic field at x, z + + Notes + ----- + As there is not plasma this is equivalent to dBz_dx + """ + return self.coilset.dB_d(x, z) + + def dBz_dx(self, x: npt.ArrayLike, z: npt.ArrayLike): + """ + Total differential of the vertical magnetic field at point (x, z) from coils + + Parameters + ---------- + x: + Radial coordinates for which to return Bz. If None, returns values + at all grid points + z: + Vertical coordinates for which to return Bz. If None, returns values + at all grid points + + Returns + ------- + : + Differential of the vertical magnetic field at x, z + + Notes + ----- + As there is not plasma this is equivalent to dBx_dz + """ + return self.coilset.dB_d(x, z) + def psi( self, x: npt.ArrayLike | None = None, @@ -905,7 +954,7 @@ class QpsiCalcMode(Enum): ZEROS = 2 -class Equilibrium(CoilSetMHDState): +class Equilibrium(CoilSetMHDState): # noqa: PLR0904 """ Represents the equilibrium state, including plasma and coil currents @@ -1557,6 +1606,46 @@ def pressure_map(self) -> npt.NDArray[np.float64]: p = self.pressure(np.clip(self.psi_norm(), 0, 1)) return p * mask + def dBx_dz(self, x: npt.ArrayLike, z: npt.ArrayLike): + """ + Total differential of the radial magnetic field at point (x, z) from coils + + Parameters + ---------- + x: + Radial coordinates for which to return Bz. If None, returns values + at all grid points + z: + Vertical coordinates for which to return Bz. If None, returns values + at all grid points + + Returns + ------- + : + Differential of the radial magnetic field at x, z + """ + return self.plasma.dBx(x, z) + self.coilset.dB_d(x, z) + + def dBz_dx(self, x: npt.ArrayLike, z: npt.ArrayLike): + """ + Total differential of the vertical magnetic field at point (x, z) from coils + + Parameters + ---------- + x: + Radial coordinates for which to return Bz. If None, returns values + at all grid points + z: + Vertical coordinates for which to return Bz. If None, returns values + at all grid points + + Returns + ------- + : + Differential of the vertical magnetic field at x, z + """ + return self.plasma.dBz(x, z) + self.coilset.dB_d(x, z) + def _get_core_mask(self) -> npt.NDArray[np.float64]: """ Get a 2-D masking array for the plasma core. diff --git a/bluemira/equilibria/plasma.py b/bluemira/equilibria/plasma.py index d9c09d60ae..20e4c51881 100644 --- a/bluemira/equilibria/plasma.py +++ b/bluemira/equilibria/plasma.py @@ -18,7 +18,12 @@ from bluemira.equilibria.constants import J_TOR_MIN from bluemira.equilibria.error import EquilibriaError from bluemira.equilibria.plotting import PlasmaCoilPlotter -from bluemira.magnetostatics.greens import greens_Bx, greens_Bz, greens_psi +from bluemira.magnetostatics.greens import ( + greens_Bx, + greens_Bz, + greens_dbz_dx, + greens_psi, +) from bluemira.utilities.tools import floatify if TYPE_CHECKING: @@ -102,15 +107,23 @@ def _set_funcs(self, plasma_psi: npt.NDArray[np.float64]): self._grid.x[:, 0], self._grid.z[0, :], plasma_psi ) self._plasma_Bx = self._Bx_func(self._grid.x, self._grid.z) + self._plasma_dBx = self._dBx_func(self._grid.x, self._grid.z) self._plasma_Bz = self._Bz_func(self._grid.x, self._grid.z) + self._plasma_dBz = self._dBz_func(self._grid.x, self._grid.z) self._plasma_Bp = np.hypot(self._plasma_Bx, self._plasma_Bz) def _Bx_func(self, x, z): return -self._psi_func(x, z, dy=1, grid=False) / x + def _dBx_func(self, x, z): + return -self._psi_func(x, z, dy=2, grid=False) / x**2 + def _Bz_func(self, x, z): return self._psi_func(x, z, dx=1, grid=False) / x + def _dBz_func(self, x, z): + return -self._psi_func(x, z, dx=2, grid=False) / x**2 + def _check_in_grid(self, x, z): return self._grid.point_inside(x, z) @@ -123,6 +136,11 @@ def _convolve(self, func, x, z): ------ EquilibriaError No known toroidal current distribution + + Returns + ------- + : + Mapped greens function """ if self._j_tor is None: raise EquilibriaError( @@ -162,7 +180,8 @@ def psi( Returns ------- - Poloidal magnetic flux at the points [V.s/rad] + : + Poloidal magnetic flux at the points [V.s/rad] """ if x is None and z is None: return self._plasma_psi @@ -193,7 +212,8 @@ def Bx( Returns ------- - Radial magnetic field at the points [T] + : + Radial magnetic field at the points [T] """ if x is None and z is None: return self._plasma_Bx @@ -202,6 +222,39 @@ def Bx( return self._convolve(greens_Bx, x, z) return self._Bx_func(x, z) + @treat_xz_array + def dBx( + self, + x: npt.ArrayLike | None = None, + z: npt.ArrayLike | None = None, + ) -> float | npt.NDArray[np.float64]: + """ + Radial magnetic field at x, z + + Parameters + ---------- + x: + Radial coordinates at which to calculate + z: + Vertical coordinates at which to calculate. + + Notes + ----- + If both x and z are None, defaults to the full map on the grid. + + Returns + ------- + : + Radial magnetic field at the points [T] + """ + if x is None and z is None: + return self._plasma_dBx + + if not self._check_in_grid(x, z): + # greens_dbx_dz not implemented, for vacuum calculations this is equivalent + return self._convolve(greens_dbz_dx, x, z) + return self._dBx_func(x, z) + @treat_xz_array def Bz( self, @@ -224,7 +277,8 @@ def Bz( Returns ------- - Vertical magnetic field at the points [T] + : + Vertical magnetic field at the points [T] """ if x is None and z is None: return self._plasma_Bz @@ -233,6 +287,38 @@ def Bz( return self._convolve(greens_Bz, x, z) return self._Bz_func(x, z) + @treat_xz_array + def dBz( + self, + x: npt.ArrayLike | None = None, + z: npt.ArrayLike | None = None, + ) -> float | npt.NDArray[np.float64]: + """ + Vertical magnetic field at x, z + + Parameters + ---------- + x: + Vertical coordinates at which to calculate + z: + Vertical coordinates at which to calculate. + + Notes + ----- + If both x and z are None, defaults to the full map on the grid. + + Returns + ------- + : + Vertical magnetic field at the points [T] + """ + if x is None and z is None: + return self._plasma_dBz + + if not self._check_in_grid(x, z): + return self._convolve(greens_dbz_dx, x, z) + return self._dBz_func(x, z) + def Bp( self, x: npt.ArrayLike | None = None, @@ -254,13 +340,14 @@ def Bp( Returns ------- - Poloidal magnetic field at the points [T] + : + Poloidal magnetic field at the points [T] """ if x is None and z is None: return self._plasma_Bp return np.hypot(self.Bx(x, z), self.Bz(x, z)) - def plot(self, ax=None): + def plot(self, ax=None) -> PlasmaCoilPlotter: """ Plot the PlasmaCoil. @@ -268,12 +355,20 @@ def plot(self, ax=None): ---------- ax: The matplotlib axes on which to plot the PlasmaCoil + + Returns + ------- + : + the plasma coil plotter object """ return PlasmaCoilPlotter(self, ax=ax) - def __repr__(self): + def __repr__(self) -> str: """ - Get a simple string representation of the PlasmaCoil. + Returns + ------- + : + A simple string representation of the PlasmaCoil. """ n_filaments = len(np.nonzero(self._j_tor > 0)[0]) return f"{self.__class__.__name__}: {n_filaments} filaments" @@ -314,7 +409,8 @@ def psi( Returns ------- - Poloidal magnetic flux at the points [V.s/rad] + : + Poloidal magnetic flux at the points [V.s/rad] """ return self._return_zeros(x, z) @@ -339,7 +435,8 @@ def Bx( Returns ------- - Radial magnetic field at the points [T] + : + Radial magnetic field at the points [T] """ return self._return_zeros(x, z) @@ -364,7 +461,8 @@ def Bz( Returns ------- - Vertical magnetic field at the points [T] + : + Vertical magnetic field at the points [T] """ return self._return_zeros(x, z) @@ -389,7 +487,60 @@ def Bp( Returns ------- - Poloidal magnetic field at the points [T] + : + Poloidal magnetic field at the points [T] + """ + return self._return_zeros(x, z) + + def dBx( + self, + x: npt.ArrayLike | None = None, + z: npt.ArrayLike | None = None, + ) -> float | npt.NDArray[np.float64]: + """ + Radial magnetic field at x, z + + Parameters + ---------- + x: + Radial coordinates at which to calculate + z: + Vertical coordinates at which to calculate. + + Notes + ----- + If both x and z are None, defaults to the full map on the grid. + + Returns + ------- + : + Radial magnetic field at the points [T] + """ + return self._return_zeros(x, z) + + def dBz( + self, + x: npt.ArrayLike | None = None, + z: npt.ArrayLike | None = None, + ) -> float | npt.NDArray[np.float64]: + """ + Vertical magnetic field at x, z + + Parameters + ---------- + x: + Vertical coordinates at which to calculate + z: + Vertical coordinates at which to calculate. + + Notes + ----- + If both x and z are None, defaults to the full map on the grid. + + Returns + ------- + : + Vertical magnetic field at the points [T] """ return self._return_zeros(x, z) diff --git a/bluemira/magnetostatics/greens.py b/bluemira/magnetostatics/greens.py index 801416a00c..598842abbf 100644 --- a/bluemira/magnetostatics/greens.py +++ b/bluemira/magnetostatics/greens.py @@ -18,6 +18,7 @@ "greens_Bx", "greens_Bz", "greens_all", + "greens_dbz_dx", "greens_dpsi_dx", "greens_dpsi_dz", "greens_psi", @@ -87,7 +88,96 @@ def ellipk_nb(k: float | np.ndarray) -> float | np.ndarray: @nb.jit(nopython=True) -def circular_coil_inductance_elliptic(radius: float, rc: float) -> float: +def _elliptic_derivatives(e, k, k2): + r"""Get :math:`\frac{dK}{dk}` and :math:`\frac{dE}{dk}` [dimensionless] + + .. math:: + + \frac{dK}{dk} &= \frac{E}{k}\frac{1}{1-k^2} - \frac{K}{k} + + \frac{dE}{dk} &= \frac{E}{k} - \frac{K}{k} + """ + sqk2 = np.sqrt(k2) + e_sqk = e / sqk2 + k_sqk = k / sqk2 + dKdk = e_sqk / (1 - k2) - k_sqk # noqa: N806 + dEdk = e_sqk - k_sqk # noqa: N806 + return dKdk, dEdk + + +@nb.jit(nopython=True) +def _dkdr(g3, xc, x): + r"""Get dkdr + + .. math:: + + \text{dkdr} &= \frac{-2 x xc \frac{x+xc}{g_3^2} + + \frac{x}{g_3}}{\sqrt{\frac{x xc}{g_3}}} + + &= -2\sqrt{x xc} \sqrt{g_3}(\frac{x+xc}{g_3^2}) + + \sqrt{\frac{g_3}{x xc}}\frac{x}{g_3} + + &= -2\sqrt{\frac{x xc}{g_3^3}} (x+xc) + \sqrt{\frac{xc}{g_3x}} + + &= -2 \frac{(x+xc)\sqrt{x xc}}{g_3^{\frac{3}{2}}} + \sqrt{\frac{xc}{g_3x}} + + unit: [m^(-1)] + """ + # old_expression = (-2 * x * xc * (x + xc) / (g3**2) + x / g3) / sqrt(x * xc / g3) + term_1 = -2 * (x + xc) * np.sqrt(x * xc) / g3**1.5 + term_2 = np.sqrt(xc / (g3 * x)) + return term_1 + term_2 + + +@nb.jit(nopython=True) +def _g(xc, zc, x, z): + r"""Get the tuple of (:math:`g_1, g_2, g_3, g_4`) + + unit: [m^2] + + .. math:: + + g_1 &= xc^2 - x^2 - z^2 + g_2 &= (xc - x)^2 + z^2 + g_3 &= (xc + x)^2 + z^2 + g_4 &= xc^2 + x^2 + z^2 + """ + x2 = x**2 + xc2 = xc**2 + z2 = (zc - z) ** 2 + g1 = xc2 - x2 - z2 + g2 = (xc - x) ** 2 + z2 + g3 = (xc + x) ** 2 + z2 + g4 = xc2 + x2 + z2 + return g1, g2, g3, g4 + + +@nb.jit(nopython=True) +def _g_r(xc, x): + r"""Get the tuple of (:math:`g_{1r}, g_{2r}, g_{3r}`) + + unit: [m] + + .. math:: + + g_{1r} &= \frac{dg_1}{dxc} = -2xc + g_{2r} &= \frac{dg_2}{dxc} = 2xc - 2x + g_{3r} &= \frac{dg_3}{dxc} = 2xc + 2x + + :math:`g_{4r}` is not used anywhere so is not computed. + """ + xc2 = 2 * xc + x2 = 2 * x + g1r = -xc2 + g2r = xc2 - x2 + g3r = xc2 + x2 + return g1r, g2r, g3r + + +@nb.jit(nopython=True) +def circular_coil_inductance_elliptic( + radius: float | np.ndarray, rc: float | np.ndarray +) -> float | np.ndarray: """ Calculate the inductance of a circular coil by elliptic integrals. @@ -255,6 +345,71 @@ def greens_dpsi_dx( return MU_0_2PI * x * ((xc**2 - (z - zc) ** 2 - x**2) * i2 + i1) +@nb.jit(nopython=True) +def greens_dbz_dx( + xc: float | np.ndarray, + zc: float | np.ndarray, + x: float | np.ndarray, + z: float | np.ndarray, +) -> float | np.ndarray: + r""" + Calculate :math:`\frac{dB_z}{dx}` (= :math:`\frac{dB_x}{dz}` for vacuum) + + Get the radial gradient of the vertical magnetic field due to a circular filament. + + unit: [N/A/m^2] + + .. math:: + + \frac{dB_z}{dx} = \frac{dB_x}{dz} = \frac{\mu_0 I}{2 \pi}\left( + - \frac{(K+E\frac{g_1}{g_2}) g_{3r}}{2 g_3^{\frac{3}{2}}} + + \frac{- E \frac{g_{1r}}{g_2} - E \frac{g_1 g_{2r}}{g_2^2} + + \frac{dE}{dk} \frac{g_1}{g_2} \text{dkdr} + + \frac{dK}{dk} \text{dkdr} + }{\sqrt{g_3}} + \right) + + where :math:`K = K(k^2), E = E(k^2)`. + + Returns + ------- + : + the gradient to the mangetic field + """ + _, k2 = calc_a_k2(xc, zc, x, z) + e, k = calc_e_k(k2) + kdk, edk = _elliptic_derivatives(e, k, k2) + g1, g2, g3, _ = _g(xc, zc, x, z) + g1r, g2r, g3r = _g_r(x, xc) + dkdr = _dkdr(g3, xc, x) + + # Avoid divide by 0 + g2 = np.where(g2 == 0, GREENS_ZERO, g2) + g3 = np.where(g3 == 0, GREENS_ZERO, g3) + logic_or = np.logical_or(g2 == 0, g3 == 0) + inv_g2 = g2**-1 + inv_g2_2 = g2**-2 + + p1 = np.where( + logic_or, + 0, + -MU_0_4PI * (k + e * g1 * inv_g2) * g3r * g3**-1.5, + ) + p2 = np.where( + logic_or, + 0, + MU_0_2PI + * ( + e * g1r * inv_g2 + - e * g1 * g2r * inv_g2_2 + + (g1 * edk * dkdr) * inv_g2 + + kdk * dkdr + ) + * g3**-0.5, + ) + return p1 + p2 + + @nb.jit(nopython=True) def greens_dpsi_dz( xc: float | np.ndarray, @@ -478,8 +633,9 @@ def greens_all( i1 *= 4 i2 *= 4 a_part = (z - zc) ** 2 + x**2 + xc**2 - b_part = -2 * x * xc - g_bx = MU_0_4PI * xc * (z - zc) * (i1 - i2 * a_part) / b_part - g_bz = MU_0_4PI * xc * ((xc + x * a_part / b_part) * i2 - i1 * x / b_part) + inv_b_part = 1 / (-2 * x * xc) + x_b_part = x * inv_b_part + g_bx = MU_0_4PI * xc * (z - zc) * (i1 - i2 * a_part) * inv_b_part + g_bz = MU_0_4PI * xc * ((xc + a_part * x_b_part) * i2 - i1 * x_b_part) g_psi = MU_0_4PI * a * ((2 - k2) * k - 2 * e) return g_psi, g_bx, g_bz From 15249ce3355472b91c0d3d8e8a39a5f7e304325e Mon Sep 17 00:00:00 2001 From: james Date: Fri, 13 Sep 2024 08:06:17 +0100 Subject: [PATCH 3/5] =?UTF-8?q?=F0=9F=8E=A8=20Coilset=20and=20equilibria?= =?UTF-8?q?=20small=20fixes?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/equilibria/coils/_grouping.py | 46 +++++++++++++++++--------- bluemira/equilibria/equilibrium.py | 1 - bluemira/equilibria/grid.py | 2 ++ bluemira/equilibria/plotting.py | 3 +- 4 files changed, 35 insertions(+), 17 deletions(-) diff --git a/bluemira/equilibria/coils/_grouping.py b/bluemira/equilibria/coils/_grouping.py index cc44ace2be..217cd2b21d 100644 --- a/bluemira/equilibria/coils/_grouping.py +++ b/bluemira/equilibria/coils/_grouping.py @@ -114,7 +114,7 @@ def __init__( Bx_analytic: bool = True, Bz_analytic: bool = True, ): - if any(not isinstance(c, Coil | CoilGroup) for c in coils): + if any(not isinstance(c, Coil | CoilGroup) for c in coils) or not coils: raise TypeError("Not all arguments are a Coil or CoilGroup.") self._coils = coils self._pad_discretisation(self.__list_getter("_quad_x")) @@ -334,7 +334,7 @@ def from_group_vecs(cls, eqdsk: EQDSKInterface) -> CoilGroup: pfcoils.append(coil) coils = pfcoils bluemira_warn( - "EQDSK coilset empty - dummy coilset in use." + "EQDSK coilset empty - dummy coilset in use. " "Please replace with an appropriate coilset." ) return cls(*coils) @@ -548,14 +548,14 @@ def _find_coil(self, name): raise KeyError(f"Coil '{name}' not found in Group") - def _get_coiltype(self, ctype: CoilType | str) -> list[Coil]: + def _get_coiltype(self, *ctype: CoilType | str) -> list[Coil]: """Find coil by type""" coils = [] - ctype = CoilType(ctype) + ctype = tuple(CoilType(ct) for ct in ctype) for c in self._coils: if isinstance(c, CoilGroup): - coils.extend(c._get_coiltype(ctype)) - elif c.ctype == ctype: + coils.extend(c._get_coiltype(*ctype)) + elif c.ctype in ctype: coils.append(c) return coils @@ -563,9 +563,17 @@ def all_coils(self) -> list[Coil]: """Get all coils as a flattened list (no CoilGroups)""" return [self[n] for n in self.name] - def get_coiltype(self, ctype: str | CoilType) -> CoilGroup | None: + def _get_type_index(self, *ctype: CoilType | str) -> npt.NDArray[int]: + coil_type = tuple(CoilType(ct) for ct in ctype) if ctype else CoilType + + return np.asarray( + [no for no, coil in enumerate(self.all_coils()) if coil.ctype in coil_type], + dtype=int, + ) + + def get_coiltype(self, *ctype: str | CoilType) -> CoilGroup | None: """Get coils matching coil type""" - if coiltype := self._get_coiltype(ctype): + if coiltype := self._get_coiltype(*ctype): return CoilGroup(*coiltype) return None @@ -1176,6 +1184,18 @@ def get_control_coils(self) -> CoilSet: coils.append(c) return CoilSet(*coils) + def get_uncontrolled_coils(self) -> CoilSet: + """Get uncontrolled coils""" + coils = [] + for c in self._coils: + if isinstance(c, CoilSet): + coils.extend(c.get_uncontrolled_coils()._coils) + elif (isinstance(c, Coil) and c.name not in self.control) or ( + isinstance(c, CoilGroup) and not any(n in self.control for n in c.name) + ): + coils.append(c) + return CoilSet(*coils) + @property def area(self) -> float: """ @@ -1213,9 +1233,9 @@ def _sum( return np.sum(output[..., inds], axis=-1) if sum_coils else output[..., inds] - def get_coiltype(self, ctype: str | CoilType) -> CoilSet | None: + def get_coiltype(self, *ctype: str | CoilType) -> CoilSet | None: """Get coils by coils type""" - if coiltype := self._get_coiltype(ctype): + if coiltype := self._get_coiltype(*ctype): return CoilSet(*coiltype) return None @@ -1229,11 +1249,7 @@ def from_group_vecs( """ self = super().from_group_vecs(eqdsk) - self.control = [ - coil.name - for ctype in control_coiltypes - for coil in self._get_coiltype(ctype) - ] + self.control = [coil.name for coil in self._get_coiltype(*control_coiltypes)] return self def get_optimisation_state( diff --git a/bluemira/equilibria/equilibrium.py b/bluemira/equilibria/equilibrium.py index 55f84fe26d..fad399f778 100644 --- a/bluemira/equilibria/equilibrium.py +++ b/bluemira/equilibria/equilibrium.py @@ -1367,7 +1367,6 @@ def solve_li( # Speed optimisations o_points, x_points = self.get_OX_points(psi=psi, force_update=True) mask = in_plasma(self.x, self.z, psi, o_points=o_points, x_points=x_points) - print() # flusher # noqa: T201 def minimise_dli(x): """ diff --git a/bluemira/equilibria/grid.py b/bluemira/equilibria/grid.py index d1a1c8c950..c0ca5d2003 100644 --- a/bluemira/equilibria/grid.py +++ b/bluemira/equilibria/grid.py @@ -58,6 +58,7 @@ class Grid: "edges", "nx", "nz", + "step", "x", "x_1d", "x_max", @@ -131,6 +132,7 @@ def __init__( [(0, z) for z in range(nz)], [(nx - 1, z) for z in range(nz)], ]) + self.step = abs(self.x_1d[0] - self.x_1d[1]) * abs(self.z_1d[0] - self.z_1d[1]) @classmethod def from_eqdict(cls, e) -> Grid: diff --git a/bluemira/equilibria/plotting.py b/bluemira/equilibria/plotting.py index 46fcb119f1..f6d2ff71d9 100644 --- a/bluemira/equilibria/plotting.py +++ b/bluemira/equilibria/plotting.py @@ -91,6 +91,7 @@ "CS": "#003688", "Plasma": "r", "NONE": "grey", + "DUM": "green", }, "edgecolor": "k", "linewidth": 1, @@ -310,7 +311,7 @@ def _get_centre(self): x, z = self._cg.get_control_coils().position xc = (max(x) + min(x)) / 2 zc = (max(z) + min(z)) / 2 - except AttributeError: + except (AttributeError, TypeError): # Not a coilset return None else: From 7f8cd514ddeed7ab1d11344f88133cbbd838d497 Mon Sep 17 00:00:00 2001 From: james Date: Fri, 13 Sep 2024 08:12:29 +0100 Subject: [PATCH 4/5] =?UTF-8?q?=E2=9C=A8=20Inductance=20additions?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/equilibria/coils/_tools.py | 43 +++++++++++++++++++---------- bluemira/magnetostatics/greens.py | 32 ++++++++++++++++++++- 2 files changed, 60 insertions(+), 15 deletions(-) diff --git a/bluemira/equilibria/coils/_tools.py b/bluemira/equilibria/coils/_tools.py index 65cbe9da10..078b09dab1 100644 --- a/bluemira/equilibria/coils/_tools.py +++ b/bluemira/equilibria/coils/_tools.py @@ -14,13 +14,19 @@ import numpy as np -from bluemira.magnetostatics.greens import circular_coil_inductance_elliptic, greens_psi +from bluemira.magnetostatics.greens import ( + circular_coil_inductance_elliptic, + greens_psi, + square_coil_inductance_kirchhoff, +) if TYPE_CHECKING: from bluemira.equilibria.coils import CoilSet -def make_mutual_inductance_matrix(coilset: CoilSet) -> np.ndarray: +def make_mutual_inductance_matrix( + coilset: CoilSet, *, square_coil: bool = False, with_quadratures: bool = False +) -> np.ndarray: """ Calculate the mutual inductance matrix of a coilset. @@ -37,28 +43,37 @@ def make_mutual_inductance_matrix(coilset: CoilSet) -> np.ndarray: ----- Single-filament coil formulation; serves as a useful approximation. """ - n_coils = coilset.n_coils() + if with_quadratures: + n_coils = coilset._quad_dx.ravel().size + xcoord = coilset._quad_x + zcoord = coilset._quad_z + dx = coilset._quad_dx + dz = coilset._quad_dz + else: + n_coils = coilset.n_coils() + xcoord = coilset.x + zcoord = coilset.z + dx = coilset.dx + dz = coilset.dz + M = np.zeros((n_coils, n_coils)) # noqa: N806 - xcoord = coilset.x - zcoord = coilset.z - dx = coilset.dx - dz = coilset.dz n_turns = coilset.n_turns itri, jtri = np.triu_indices(n_coils, k=1) - M[itri, jtri] = ( n_turns[itri] * n_turns[jtri] - * greens_psi(xcoord[itri], zcoord[itri], xcoord[jtri], zcoord[jtri]) + * greens_psi(xcoord[itri], zcoord[itri], xcoord[jtri], zcoord[jtri]).ravel() ) M[jtri, itri] = M[itri, jtri] - radius = np.hypot(dx, dz) - for i in range(n_coils): - M[i, i] = n_turns[i] ** 2 * circular_coil_inductance_elliptic( - xcoord[i], radius[i] - ) + ind = np.diag_indices_from(M) + + if square_coil: + M[ind] = square_coil_inductance_kirchhoff(xcoord, dx, dz) + else: + radius = np.hypot(dx, dz) + M[ind] = n_turns**2 * circular_coil_inductance_elliptic(xcoord, radius) return M diff --git a/bluemira/magnetostatics/greens.py b/bluemira/magnetostatics/greens.py index 598842abbf..b23fa7ec66 100644 --- a/bluemira/magnetostatics/greens.py +++ b/bluemira/magnetostatics/greens.py @@ -215,7 +215,9 @@ def circular_coil_inductance_elliptic( return MU_0 * (2 * radius - rc) * ((1 - k**2 / 2) * ellipk_nb(k) - ellipe_nb(k)) -def circular_coil_inductance_kirchhoff(radius: float, rc: float) -> float: +def circular_coil_inductance_kirchhoff( + radius: float | np.ndarray, rc: float | np.ndarray +) -> float | np.ndarray: """ Calculate the inductance of a circular coil by Kirchhoff's approximation. @@ -239,6 +241,34 @@ def circular_coil_inductance_kirchhoff(radius: float, rc: float) -> float: return MU_0 * radius * (np.log(8 * radius / rc) - 2 + 0.25) +def square_coil_inductance_kirchhoff( + radius: float | np.ndarray, width: float | np.ndarray, height: float | np.ndarray +) -> float | np.ndarray: + """ + Calculate the inductance of a square coil by Kirchhoff's approximation. + + radius: + The radius of the square coil + width: + The width of the coil cross-section + height + The height of the coil cross-section + + Returns + ------- + The self-inductance of the square coil [H] + + Notes + ----- + .. math:: + + Inductance = \\mu_0 radius (ln(8\\frac{radius}{width + height}) - 0.5) + + where :math:`\\mu_{0}` is the vacuum permeability + """ + return MU_0 * radius * (np.log(8 * radius / (width + height)) - 0.5) + + @nb.jit(nopython=True) def greens_psi( xc: float | np.ndarray, From df2e52643bd4c56ddcaab063975597e872aad712 Mon Sep 17 00:00:00 2001 From: james Date: Fri, 13 Sep 2024 08:13:36 +0100 Subject: [PATCH 5/5] =?UTF-8?q?=E2=9C=A8=20Rzip?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bluemira/equilibria/vertical_stability.py | 302 ++++++++++++++++++++++ 1 file changed, 302 insertions(+) create mode 100644 bluemira/equilibria/vertical_stability.py diff --git a/bluemira/equilibria/vertical_stability.py b/bluemira/equilibria/vertical_stability.py new file mode 100644 index 0000000000..551983100f --- /dev/null +++ b/bluemira/equilibria/vertical_stability.py @@ -0,0 +1,302 @@ +# 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 + +""" +Vertical stability calculations +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +import numpy as np + +from bluemira.base.look_and_feel import bluemira_debug, bluemira_warn +from bluemira.equilibria.coils._tools import make_mutual_inductance_matrix + +if TYPE_CHECKING: + import numpy.typing as npt + + from bluemira.equilibria.coils._grouping import CoilSet + from bluemira.equilibria.equilibrium import Equilibrium + + +def calculate_rzip_stability_criterion(eq: Equilibrium) -> float: + """ + Calculate the rzip stability criterion for a given equilibrium and coilset + + Returns + ------- + : + Stability criterion + + Notes + ----- + See ~:class:`~bluemira.equilibria.vertical_stability.RZIp` for details + """ + return RZIp(eq.coilset)(eq) + + +class RZIp: + """ + RZIp model + + Parameters + ---------- + coilset: + The full coilset + + Notes + ----- + A value ~ 1.5 considered optimal and controllable. + See https://doi.org/10.13182/FST89-A39747 for further explanation + + < 1 plasma mass becomes a factor + = 1 massless plasma solution not valid (said to represent MHD effects) + > 1 displacement growth dominated by L/R of passive system + + .. math:: + + f = -\\frac{F_s}{F_d} + = \\frac{I_p^T M^{\\prime}_p|_s [M_s|_s]^{-1} M^{\\prime}_s|_p I_p} + {I_p^T M^{\\prime\\prime}_p|_c I_c} + """ + + def __init__(self, coilset): + self.coilset = coilset + + @property + def coilset(self): + return self._coilset + + @coilset.setter + def coilset(self, cs: CoilSet): + # TODO @je-cook Possibly implement quad indexing? + self._coilset = cs + self.ind_mat = make_mutual_inductance_matrix(cs, square_coil=True) + + def __call__( + self, + eq: Equilibrium, + cc_current_vector: None | npt.NDArray = None, + ) -> float: + if not hasattr(eq, "profiles"): + bluemira_warn("No profiles found on equilibrium") + # this constraint is pretty irrelevant for breakdown + return 0 + + cc = self.coilset.get_control_coils() + + if cc_current_vector is None: + cc_current_vector = cc.current + else: + cc.current = cc_current_vector + + eq._remap_greens() + + control_ind = self.coilset._control_ind + uncontrolled_ind = list( + set(self.coilset._get_type_index()) - set(self.coilset._control_ind) + ) + + return stab_destab( + cc_current=cc_current_vector, + # res_ps=self.coilset.get_uncontrolled_coils().resistance, + l_ps_ps=self.ind_mat[uncontrolled_ind][:, uncontrolled_ind], + l_ac_ps=self.ind_mat[control_ind][:, uncontrolled_ind], + l_ac_ac=self.ind_mat[control_ind][:, control_ind], + # ind=M, + r_struct=self.coilset.x, + x_ac=cc.x, + i_plasma=eq._jtor * eq.grid.step, + br_struct_grid=eq._bx_green, + dbrdz_struct_grid=eq._db_green[..., control_ind], + # max_eigens=self.max_eigens, + ) + + +def stab_destab( + cc_current, + # res_ps, + l_ps_ps, + l_ac_ps, + l_ac_ac, + r_struct, + x_ac, + i_plasma, + br_struct_grid, + dbrdz_struct_grid, + # max_eigens, +) -> float: + """ + Calculate the stabilising / destabilising effect of the equilibria and structures + """ + # # single filament approx l_acps_eigen_vec + # eigen_values, eigen_vectors = eigenv(res_ps, l_ps_ps, max_eigens=max_eigens) + # l_acps_eigen_vec = l_ac_ps @ eigen_vectors + + # # l_struct in rzip + # mss = np.hstack([ + # np.vstack([l_ac_ac, l_acps_eigen_vec.T]), + # np.vstack([l_acps_eigen_vec, eigen_values]), + # ]) + + # l_struct in rzip + mss = np.vstack([ + np.hstack([l_ac_ac, l_ac_ps]), + np.hstack([l_ac_ps.T, l_ps_ps]), + ]) + + msp_prime = (2 * np.pi * r_struct * br_struct_grid).reshape(-1, r_struct.size) + i_plasma = i_plasma.reshape(-1) + + # stabilising force/ destabilising force differentiated wrt to z coord + # f = -d_fs / d_fd + destabilising = np.einsum( + "b, bd, d", + i_plasma, + (2 * np.pi * x_ac * dbrdz_struct_grid).reshape(-1, cc_current.size), + cc_current, + optimize=["einsum_path", (0, 1), (0, 1)], + ) + + stabilising = np.einsum( + "d, db, bc, ac, a", + i_plasma, + msp_prime, + np.linalg.inv(mss), + msp_prime, + i_plasma, + optimize=["einsum_path", (0, 1), (1, 2), (0, 1), (0, 1)], + ) + + # from scipy.io import loadmat + + # ml = loadmat("leuer_test.mat") + # import ipdb + + # ipdb.set_trace() + # not infinite if destabilising is 0 because therefore it is stable + criterion = (-stabilising / destabilising) if destabilising != 0 else 0 + bluemira_debug(f"{stabilising=}, {destabilising=}, {criterion=}") + return criterion + + +# def mutual_ind(xc, zc, x, z): +# r"""File to compute the mutual inductance between two coaxial circular filaments, +# i.e. the planes of the two coils are parallel and they share the same central axis. +# (William R. Smythe: Static and Dynaimc Electricity Third Edition, p.335 equation 1) + +# .. math:: + +# M_{12} &= 2\mu_0 k^{-1} (ab)^{1/2} \left[(1-\frac{1}{2}k^2)K-E\right] + +# &= \frac{2 \mu_0 }{k} \sqrt{ab}\left[\left(1-\frac{k^2}{2}\right)K-E\right] + +# where + +# .. math:: + +# k^2 = \frac{4ab}{(a+b)^2 + z^2} + +# Note the symmetry between :math:`a` and :math:`b` in the equations: +# :math:`a` and :math:`b` are commutable. + +# Parameters +# ---------- +# a, b : radii of coil 1 and coil 2 respectively [m] +# (denoted with :math:`a` and :math:`b` in the equation) + +# z: separation between the planes of coils [m] + +# Returns +# ------- +# Mutual Inductance: unit: [N/m/A^2] +# """ +# _, k2 = calc_a_k2(xc, zc, x, z) +# e, k = calc_e_k(k2) +# i_ch = np.nonzero(k2 == 1) +# calc_a_k2(xc, zc, x, z) +# # TODO: review this hack! +# k2[i_ch] = 0.1 +# mut_ind = 2 * MU_0 / np.sqrt(k2) * np.sqrt(xc * zc) * ((1 - 0.5 * k2) * k - e) +# mut_ind[i_ch] = 0 +# return mut_ind + + +def eigenv(res_ps, l_ps_ps, max_eigens=40): + inv_r = np.linalg.inv(res_ps[..., None] * np.eye(res_ps.shape[-1])) + inv_r_l = np.einsum("...ab, ...bc -> ...ac", inv_r, l_ps_ps) + eig_val, eig_vec = np.linalg.eig(np.asarray(inv_r_l, dtype=np.complex128)) + + # sort (eig_vec_sort) into ascending importance of values (vectors Es and values Vs) + sort_ind = np.argsort(eig_val, axis=-1)[..., ::-1] + eig_val = eig_val[..., sort_ind].real + eig_vec = eig_vec[..., sort_ind].real + norm_res = 1 / np.sum( # TODO compare with rzip with multiple structures + 1 / np.atleast_2d(res_ps), axis=-1 + ) + norm_mat = np.einsum( # diagonalise + "...aa -> ...a", + np.einsum( # arr1 @ arr2 @ inv(arr3.T) + "...ab, ...bc, ...cd -> ...ad", + norm_res * np.linalg.inv(eig_vec), + inv_r, + np.linalg.inv(np.einsum("...ij -> ...ji", eig_vec)), + optimize=["einsum_path", (0, 1), (0, 1)], + ), + ) + idx = np.nonzero(norm_mat >= 0) + norm_mat[idx] = np.sqrt(norm_mat[idx]) + idx = np.nonzero(norm_mat < 0) + norm_mat[idx] = -np.sqrt(-norm_mat[idx]) + + eig_val *= norm_res + eig_vec @= norm_mat[..., None] * np.eye(norm_mat.shape[-1]) + + # eig_val = eig_val[..., :max_eigens] + # eig_vec = eig_vec[..., :max_eigens] + return eig_val * np.eye(eig_val.shape[0]), eig_vec + + # TODO inter structure interactions + eigen_vectors_2 = np.zeros((n_passive_fils, n_eigen_modes_max)) + eigen_values_2 = np.zeros((1, n_eigen_modes_max)) + + for k in range(n_passive_struct): + passive_range = slice(passive_coil_ranges[k, 0], passive_coil_ranges[k, 1]) + eigen_range_start = n_eigen_modes_array_max[k] if k > 0 else 0 + eigen_range = slice( + eigen_range_start, n_eigen_modes_array_max[k] + eigen_range_start + ) + + # structure to structure terms + # TODO untested + for k2 in range(k + 1, n_passive_struct): + raise NotImplementedError(">1 passive structure not implemented") + passive_range2 = slice( + passive_coil_ranges[k2, 0], passive_coil_ranges[k2, 1] + ) + + eigen_range_start2 = sum(n_eigen_modes_array_max[:k2]) + eigen_range2 = slice( + eigen_range_start2, n_eigen_modes_array_max[k2] + eigen_range_start2 + ) + + et_me = ( + eigen_vectors[k2, :, : n_eigen_modes_array_max[k]] + @ l_pf_pf[passive_range, passive_range2] + ) @ eigen_vectors[k2, :, : n_eigen_modes_array_max[k2]] + eigen_values_2[k, eigen_range, eigen_range2] = et_me + eigen_values_2[k, eigen_range2, eigen_range] = et_me + + # TODO if there is more than 1 passive structure it will break + eigen_values_2[eigen_range] = eigen_values[k] + eigen_vectors_2[passive_range, eigen_range] = eigen_vectors[k] + + eigen_vectors = eigen_vectors_2 + eigen_values = eigen_values_2 + + return eigen_values, eigen_vectors