Add VerticalFieldConstraint and function to estimate plasma equilibrium field

bluemira/equilibria/optimisation/constraints.py

@@ -550,6 +550,49 @@
         return 2
 
 
+class VerticalFieldConstraint(AbsoluteMagneticConstraint):
+    """
+    Absolute vertical field (Bz) constraint.
+    """
+
+    def __init__(
+        self,
+        x: Union[float, np.ndarray],
+        z: Union[float, np.ndarray],
+        target_value: float,
+        weights: Union[float, np.ndarray] = 1.0,
+        tolerance: Union[float, np.ndarray] = 1e-6,
+    ):
+        super().__init__(
+            x,
+            z,
+            target_value,
+            weights=weights,
+            tolerance=tolerance,
+            f_constraint=AxBConstraint,
+            constraint_type="equality",
+        )
+
+    def control_response(self, coilset: CoilSet) -> np.ndarray:
+        """
+        Calculate control response of a CoilSet to the constraint.
+        """
+        return coilset.Bz_response(self.x, self.z, control=True)
+
+    def evaluate(self, eq: Equilibrium) -> np.ndarray:
+        """
+        Calculate the value of the constraint in an Equilibrium.
+        """
+        return eq.Bz(self.x, self.z)
+
+    def plot(self, ax):
+        """
+        Plot the constraint onto an Axes.
+        """
+        kwargs = {"marker": "^", "markersize": 8, "color": "b", "linestyle": "None"}
+        ax.plot(self.x, self.z, **kwargs)
+
+
 class PsiConstraint(AbsoluteMagneticConstraint):
     """
     Absolute psi value constraint.

bluemira/plasma_physics/rules_of_thumb.py

@@ -8,9 +8,11 @@
 A collection of simple 0-D rules of thumb for tokamak plasmas.
 """
 
+from typing import Optional
+
 import numpy as np
 
-from bluemira.base.constants import EV_TO_J, K_BOLTZMANN, MU_0
+from bluemira.base.constants import EV_TO_J, K_BOLTZMANN, MU_0, MU_0_4PI
 from bluemira.plasma_physics.collisions import coulomb_logarithm, spitzer_conductivity
 
 
@@ -282,3 +284,41 @@
     """
     nu = q_star / q_0 - 1.0
     return np.log(1.65 + 0.89 * nu)
+
+
+def estimate_vertical_field(
+    R_0: float,
+    A: float,
+    I_p: float,
+    beta_p_th: float,
+    l_i: float,
+    kappa_95: Optional[float] = None,
+) -> float:
+    """
+    Estimate the vertical field to keep the plasma in equilibrium.
+
+    Parameters
+    ----------
+    R_0:
+        Plasma major radius [m]
+    A:
+        Plasma aspect ratio
+    I_p:
+        Plasma current [A]
+    beta_p_th:
+        Thermal poloidal beta
+    l_i:
+        Normalised internal inductance
+    kappa:
+        Plasma elongation at the 95th percentile flux surface
+
+    Notes
+    -----
+    See e.g. Ferrara et al., "Alcasim simulation code for Alcator C-Mod"
+    https://ieeexplore.ieee.org/stamp/stamp.jsp?tp=&arnumber=4178094
+
+    The kappa term is not always present in textbooks and the like, and
+    is almost certainly irrelevant at the end of breakdown.
+    """
+    k_term = 1.0 if kappa_95 is None else np.sqrt(2.0 / (1.0 + kappa_95**2))
+    return -MU_0_4PI * I_p / R_0 * (beta_p_th + 0.5 * l_i - 1.5 + np.log(8 * A * k_term))

tests/equilibria/test_optimisation.py

@@ -5,17 +5,21 @@
 # SPDX-License-Identifier: LGPL-2.1-or-later
 import numpy as np
 
-from bluemira.equilibria import Equilibrium
+from bluemira.equilibria import Breakdown, Equilibrium
 from bluemira.equilibria.coils import Coil, CoilSet
 from bluemira.equilibria.grid import Grid
 from bluemira.equilibria.optimisation.constraints import (
     FieldNullConstraint,
     IsofluxConstraint,
     MagneticConstraintSet,
+    VerticalFieldConstraint,
+)
+from bluemira.equilibria.optimisation.problem import (
+    TikhonovCurrentCOP,
 )
-from bluemira.equilibria.optimisation.problem import TikhonovCurrentCOP
 from bluemira.equilibria.profiles import CustomProfile
 from bluemira.equilibria.solve import PicardIterator
+from bluemira.optimisation import Algorithm
 
 
 def coilset_setup(materials=False):
@@ -79,3 +83,24 @@ def test_isoflux_constrained_tikhonov_current_optimisation():
         ],
         decimal=5,
     )
+
+
+def test_vertical_field_constraint():
+    coilset = coilset_setup()
+    grid = Grid(4.5, 14, -9, 9, 65, 65)
+    eq = Breakdown(coilset, grid)
+
+    bz_constraint = VerticalFieldConstraint(9, 0.0, -0.75, 1.0, tolerance=1e-6)
+    targets = MagneticConstraintSet([bz_constraint])
+    opt_problem = TikhonovCurrentCOP(
+        eq.coilset,
+        eq,
+        targets,
+        gamma=1e-18,
+        opt_algorithm=Algorithm.SLSQP,
+        opt_conditions={"max_eval": 5},
+    )
+
+    coilset = opt_problem.optimise()
+
+    np.testing.assert_allclose(eq.Bz(9.0, 0.0), -0.75, rtol=1e-6)