diff --git a/examples/data/simple_diverted_currents_PaxisIp.pk b/examples/data/simple_diverted_currents_PaxisIp.pk
index eae8a38d..08a31ddf 100644
Binary files a/examples/data/simple_diverted_currents_PaxisIp.pk and b/examples/data/simple_diverted_currents_PaxisIp.pk differ
diff --git a/examples/data/simple_limited_currents_PaxisIp.pk b/examples/data/simple_limited_currents_PaxisIp.pk
index 148f88ae..c3ec9ad9 100644
Binary files a/examples/data/simple_limited_currents_PaxisIp.pk and b/examples/data/simple_limited_currents_PaxisIp.pk differ
diff --git a/examples/example01b - advanced_static_inverse_solve.ipynb b/examples/example01b - advanced_static_inverse_solve.ipynb
index 95f6b332..1a75390c 100644
--- a/examples/example01b - advanced_static_inverse_solve.ipynb
+++ b/examples/example01b - advanced_static_inverse_solve.ipynb
@@ -70,7 +70,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "### (Normalised) poloidal flux constraints at specific locations\n",
+ "### Feature 1: (Normalised) poloidal flux constraints at specific locations\n",
"\n",
"In addition to the `null points`, `isoflux set`, and `coil_current_limits` constraints introduced in the previous notebook, additional methods are available for more advanced control of the inverse problem.\n",
"\n",
@@ -223,7 +223,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "#### Weighting constraints within the solver\n",
+ "### Feature 2: Weighting constraints within the solver\n",
"\n",
"It is possible to specify **weights** on the different types of constraints within the inverse solver, enabling the solver to prioritise certain constraints over others.\n",
"\n",
@@ -411,12 +411,283 @@
"plt.tight_layout()"
]
},
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "### Feature 3: second-order magnetic null constraints (i.e. generating snowflake divertors)\n",
+ "\n",
+ "Now we'll try to use the more complex constraints to form second order nulls (using the `null_points_2nd_order` parameter) and a snowflake divertor geometry. \n",
+ "\n",
+ "\n",
+ "In some situations, it may be desirable to control the **local structure of the magnetic field** around a null point. This is achieved by enforcing **second-order null point constraints**, which additionally constrain the *spatial derivatives* of the magnetic field.\n",
+ "\n",
+ "Recall that the poloidal magnetic field is related to the poloidal flux $\\psi$ via\n",
+ "\n",
+ "$$\n",
+ "B_R = -\\frac{1}{R} \\frac{\\partial \\psi}{\\partial Z}, \n",
+ "\\qquad\n",
+ "B_Z = \\frac{1}{R} \\frac{\\partial \\psi}{\\partial R}.\n",
+ "$$\n",
+ "\n",
+ "A standard null point $(R_X, Z_X)$ satisfies\n",
+ "\n",
+ "$$\n",
+ "B_R(R_X, Z_X) = B_Z(R_X, Z_X) = 0.\n",
+ "$$\n",
+ "\n",
+ "Second-order null point constraints extend this by additionally imposing conditions on the **Jacobian of the magnetic field**, i.e. its first spatial derivatives:\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial B_R}{\\partial R}, \\quad\n",
+ "\\frac{\\partial B_R}{\\partial Z}, \\quad\n",
+ "\\frac{\\partial B_Z}{\\partial R}, \\quad\n",
+ "\\frac{\\partial B_Z}{\\partial Z}.\n",
+ "$$\n",
+ "\n",
+ "In FreeGSNKE, the following constraints are enforced at each second-order null point:\n",
+ "\n",
+ "$$\n",
+ "B_R = 0, \\qquad B_Z = 0,\n",
+ "$$\n",
+ "\n",
+ "$$\n",
+ "\\frac{\\partial B_R}{\\partial R} = 0, \\qquad\n",
+ "\\frac{\\partial B_Z}{\\partial Z} = 0, \\qquad\n",
+ "\\frac{\\partial B_R}{\\partial Z} = 0, \\qquad\n",
+ "\\frac{\\partial B_Z}{\\partial R} = 0.\n",
+ "$$\n",
+ "\n",
+ "This results in a total of **six constraints per null point** so be aware that:\n",
+ "- imposing too many (second-order) constraints can easily **overconstrain** the system.\n",
+ "- in practice, they are most useful when combined with a limited number of other constraints (e.g. isoflux sets).\n",
+ "\n",
+ "Here, we'll use one second-order null constraint in combination with a few isoflux constraints. \n",
+ "\n",
+ "Also note how we increase `constrain.mu_coils` as snowflake geometries tend to demand a lot from the PF coils. We need to make sure we don't violate their limits. \n",
+ "\n",
+ "But first lets reset the equilibrium and profile objects (cranking up the resolution slightly)."
+ ]
+ },
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
- "source": []
+ "source": [
+ "# we are modifying the control, build a new tokamak\n",
+ "tokamak = build_machine.tokamak(\n",
+ " active_coils_path=\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n",
+ " passive_coils_path=\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n",
+ " limiter_path=\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n",
+ " wall_path=\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n",
+ ")\n",
+ "\n",
+ "# a new eq object resets the coil currents and equilibrium\n",
+ "eq = equilibrium_update.Equilibrium(\n",
+ " tokamak=tokamak, # provide tokamak object\n",
+ " Rmin=0.1, Rmax=2.0, # radial range\n",
+ " Zmin=-2.2, Zmax=2.2, # vertical range\n",
+ " nx=129, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n",
+ " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n",
+ ")\n",
+ "\n",
+ "# reset the profiles\n",
+ "profiles = ConstrainPaxisIp(\n",
+ " eq=eq, # equilibrium object\n",
+ " paxis=8e3, # profile object\n",
+ " Ip=6e5, # plasma current\n",
+ " fvac=0.5, # fvac = rB_{tor}\n",
+ " alpha_m=1.8, # profile function parameter\n",
+ " alpha_n=1.2 # profile function parameter\n",
+ ")\n",
+ "\n",
+ "# re-initialise the solver object (due to new grid resolution)\n",
+ "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n",
+ "\n",
+ "# first we specify some alternative constraints\n",
+ "Rout = 1.32 # outboard midplane radius\n",
+ "Rin = 0.3 # inboard midplane radius\n",
+ "\n",
+ "# locations of X-points (this time they will be second order nulls!)\n",
+ "Rx = 0.65\n",
+ "Zx = 1.2\n",
+ "\n",
+ "# second order null constraints to form the snowflake\n",
+ "null_points_2nd_order = [[Rx], [Zx]]\n",
+ "\n",
+ "# isoflux constraints (no divertor constraints this time)\n",
+ "isoflux_set = np.array([[[Rx, Rin, Rout], [Zx, 0.,0.]]])\n",
+ "\n",
+ "# starting guess for the Solenoid current\n",
+ "# eq.tokamak.set_coil_current('Solenoid', 5000)\n",
+ "# eq.tokamak['Solenoid'].control = True\n",
+ "\n",
+ "# also make it perfectly up/down symmetric\n",
+ "eq.tokamak.set_coil_current('P6', 0)\n",
+ "eq.tokamak['P6'].control = False\n",
+ "\n",
+ "coil_current_limits = [\n",
+ "# upper limits...\n",
+ "# P1, PX, D1, D2, D3, Dp, D5, D6, D7, P4, P5\n",
+ " [1e4, 6e3, 9e3, 9e3, 7e3, 7e3, 5e3, 3.5e3, 5e3, 0, 0],\n",
+ "# lower limits...\n",
+ "# P1, PX, D1, D2, D3, Dp, D5, D6, D7, P4, P5\n",
+ " [-1e4, -6e3, -9e3, -9e3, -7e3, -7e3, -5e3, -3.5e3, -5e3, -1.2e4, -1.2e4]\n",
+ "]\n",
+ "\n",
+ "# pass the magnetic constraints to a new constrain object\n",
+ "constrain = Inverse_optimizer(\n",
+ " null_points_2nd_order=null_points_2nd_order,\n",
+ " isoflux_set=isoflux_set,\n",
+ " coil_current_limits=coil_current_limits,\n",
+ ")\n",
+ "\n",
+ "# increase this to ensure coil currents stay with limits\n",
+ "constrain.mu_coils = 5e7\n",
+ "\n",
+ "# carry out the solve\n",
+ "GSStaticSolver.solve(eq=eq, \n",
+ " profiles=profiles, \n",
+ " constrain=constrain, \n",
+ " target_relative_tolerance=1e-6,\n",
+ " target_relative_psit_update=1e-3,\n",
+ " verbose=True,\n",
+ " l2_reg=np.array([1e-12]*11),\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot the resulting equilbrium\n",
+ "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=100)\n",
+ "ax1.grid(True, which='both')\n",
+ "eq.plot(axis=ax1, show=False)\n",
+ "eq.tokamak.plot(axis=ax1, show=False)\n",
+ "constrain.plot(axis=ax1,show=True)\n",
+ "ax1.set_xlim(0.1, 2.15)\n",
+ "ax1.set_ylim(-2.25, 2.25)\n",
+ "plt.tight_layout()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Zoom in to see the snowflake divertor geometry!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# upper X-point\n",
+ "fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8), dpi=80)\n",
+ "ax1.grid(True, which='both')\n",
+ "eq.plot(axis=ax1, show=False)\n",
+ "eq.tokamak.plot(axis=ax1, show=False)\n",
+ "constrain.plot(axis=ax1,show=False)\n",
+ "ax1.set_xlim(0.15, 1.15)\n",
+ "ax1.set_ylim(0.75, 1.75)\n",
+ "plt.tight_layout()\n",
+ "\n",
+ "# lower X-point\n",
+ "fig1, ax1 = plt.subplots(1, 1, figsize=(8, 8), dpi=80)\n",
+ "ax1.grid(True, which='both')\n",
+ "eq.plot(axis=ax1, show=False)\n",
+ "eq.tokamak.plot(axis=ax1, show=False)\n",
+ "constrain.plot(axis=ax1,show=False)\n",
+ "ax1.set_xlim(0.15, 1.15)\n",
+ "ax1.set_ylim(-1.75, -0.75)\n",
+ "plt.tight_layout()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "Just double check that the coil currents stayed within the limits!"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# storage\n",
+ "coil_names = []\n",
+ "active_currents = []\n",
+ "min_limits = []\n",
+ "max_limits = []\n",
+ "\n",
+ "for coil_name, coil_current, ul, ll in zip(\n",
+ " eq.tokamak.coils_dict,\n",
+ " eq.tokamak.getCurrentsVec()[:12],\n",
+ " coil_current_limits[0] + [None], # upper limits\n",
+ " coil_current_limits[1] + [None], # lower limits\n",
+ "):\n",
+ " coil_names.append(coil_name)\n",
+ " active_currents.append(coil_current)\n",
+ " min_limits.append(ll)\n",
+ " max_limits.append(ul)\n",
+ "\n",
+ " # skip coils without limits\n",
+ " if ul is not None or ll is not None:\n",
+ " if (coil_current >= ll) & (coil_current <= ul):\n",
+ " within_limit = True\n",
+ " else:\n",
+ " within_limit = False\n",
+ " else:\n",
+ " within_limit = None\n",
+ "\n",
+ " # print\n",
+ " print(f\"{coil_name}: {coil_current:.2f} --> [{ll},{ul}] [A] --> limits met = {within_limit}\")\n",
+ "\n",
+ "# convert to numpy\n",
+ "active_currents = np.array(active_currents)\n",
+ "min_limits = np.array(min_limits)\n",
+ "max_limits = np.array(max_limits)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# VISUALISE HOW CLOSE THEY ARE TO THE LIMITS (normalised)\n",
+ "\n",
+ "# normalized to [-1, 1] scale (relative to min/max)\n",
+ "min_limits[min_limits==None] = np.inf\n",
+ "max_limits[max_limits==None] = np.inf\n",
+ "norm_current = 2 * (active_currents - min_limits) / (max_limits - min_limits) - 1\n",
+ "\n",
+ "# plot\n",
+ "plt.figure(figsize=(12, 5), dpi=70)\n",
+ "x = np.arange(len(coil_names))\n",
+ "\n",
+ "# line plot\n",
+ "plt.plot(x, norm_current, marker='x', linestyle='-', color='red', label=\"Current\")\n",
+ "plt.axhline(-1, color='k', linestyle='-', linewidth=1, label=\"Min. limit\")\n",
+ "plt.axhline(1, color='k', linestyle='-', linewidth=1, label=\"Max. limit\")\n",
+ "\n",
+ "# formatting\n",
+ "plt.xticks(x, coil_names, rotation=45, ha='right')\n",
+ "plt.ylabel('Coil name [Amps]')\n",
+ "plt.ylabel('Normalised current')\n",
+ "plt.title('Coil current proximity to limits (normalised)')\n",
+ "plt.grid(True, linestyle='--', alpha=0.6)\n",
+ "plt.legend()\n",
+ "plt.tight_layout()\n",
+ "plt.show()\n"
+ ]
},
{
"cell_type": "code",
diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py
index aaa3ff32..2866920c 100644
--- a/freegsnke/inverse.py
+++ b/freegsnke/inverse.py
@@ -14,9 +14,9 @@
it under the terms of the GNU Lesser General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
-
+
You should have received a copy of the GNU Lesser General Public License
-along with FreeGSNKE. If not, see .
+along with FreeGSNKE. If not, see .
"""
import itertools
@@ -28,13 +28,14 @@
class Inverse_optimizer:
"""This class implements a gradient based optimiser for the coil currents,
- used to perform (static) inverse GS solves.
+ used to perform (static) inverse Grad-Shafranov solves.
"""
def __init__(
self,
isoflux_set=None,
null_points=None,
+ null_points_2nd_order=None,
psi_vals=None,
coil_current_limits=None,
psi_norm_limits=None,
@@ -73,8 +74,8 @@ def __init__(
where:
Rcoords : 1D array of radial coordinates
Zcoords : 1D array of vertical coordinates
- weights : (optional, array of 1's if not provided)1D array of weights to increase/decrease
- the influence of the constraint on the solution.
+ weights : (optional, array of 1's if not provided) 1D array of weights
+ to increase/decrease the influence of the constraint on the solution.
All specified points within each set are required to share
the same poloidal flux value.
@@ -89,6 +90,17 @@ def __init__(
• X-points
• O-points (magnetic axes)
+ null_points_2nd_order : list or ndarray, optional
+ Second-order magnetic null point constraints.
+
+ Structure:
+ [Rcoords, Zcoords]
+
+ Specifies coordinates of desired 2nd-order null points, typically
+ used for snowflake divertor configurations.
+
+ Note: Do not repeat coordinates already provided in ``null_points``.
+
psi_vals : list or ndarray, optional
Direct flux value constraints.
@@ -106,14 +118,10 @@ def __init__(
Structure:
[upper_limits, lower_limits]
- Each entry is a list with length equal to the number of
- controllable coils.
+ Each entry is a list with length equal to the number of controllable coils.
Example:
- [
- [Imax1, Imax2, ...],
- [Imin1, Imin2, ...]
- ]
+ [ [Imax1, Imax2, ...], [Imin1, Imin2, ...] ]
Use None to indicate no bound.
@@ -124,33 +132,37 @@ def __init__(
[Rcoord, Zcoord, normalised_psi_value, constraint_sign]
Constraint form:
-
- If constraint_sign = 1:
- ψ_norm ≥ ψ_target
-
- If constraint_sign = -1:
- ψ_norm ≤ ψ_target
+ If constraint_sign = 1: ψ_norm ≥ ψ_target
+ If constraint_sign = -1: ψ_norm ≤ ψ_target
Normalised flux is defined:
ψ_norm = (ψ - ψ_axis) / (ψ_boundary - ψ_axis)
weight_isoflux : float
- The weight of the isoflux constraints in the least-squares optimisation problem (default = 1.0).
+ The weight of the isoflux constraints in the least-squares optimisation
+ problem (default = 1.0).
+
weight_nulls : float
- The weight of the null point (X-point) constraints in the least-squares optimisation problem (default = 1.0).
+ The weight of the null point (X-point) constraints in the least-squares
+ optimisation problem (default = 1.0).
+
weight_psi : float
- The weight of the psi value constraints in the least-squares optimisation problem (default = 1.0).
+ The weight of the psi value constraints in the least-squares optimisation
+ problem (default = 1.0).
+
mu_coils : float
- A penalty factor applied to violation of the coil current limits (default = 1e5).
+ A penalty factor applied to violation of the coil current limits
+ (default = 1e5).
+
mu_psi_norm : float
- A penalty factor applied to violation of the normalised psi limits (default = 1e5).
+ A penalty factor applied to violation of the normalised psi limits
+ (default = 1e6).
Notes
-----
- Increasing the weights/penalty factors causes the least-squares optimisation to proritise satisfying the
- higher-weighted/penaltied constraints.
+ Increasing the weights/penalty factors causes the least-squares optimisation
+ to prioritise satisfying the higher-weighted/penalised constraints.
"""
-
# ------------------------------------------------------------
# Isoflux constraint processing
# ------------------------------------------------------------
@@ -187,10 +199,18 @@ def __init__(
if self.null_points is not None:
self.null_points = np.array(self.null_points)
+ # ------------------------------------------------------------
+ # Second-order null point constraints (snowflakes)
+ # ------------------------------------------------------------
+ self.null_points_2nd_order = null_points_2nd_order
+ if self.null_points_2nd_order is not None:
+ self.null_points_2nd_order = np.array(self.null_points_2nd_order)
+
# ------------------------------------------------------------
# Direct flux value constraints
# These impose ψ(R,Z) = ψ_target at specified locations
# ------------------------------------------------------------
+
self.psi_vals = psi_vals
if self.psi_vals is not None:
@@ -563,9 +583,34 @@ def build_greens(self, eq):
)
# ------------------------------------------------------------
- # Flux value constraint Green's functions
+ # Magnetic null point Green's functions (for second-order nulls)
# ------------------------------------------------------------
+ if self.null_points_2nd_order is not None:
+ # for first order null
+ self.Gbr_2nd_order = eq.tokamak.createBrGreensVec(
+ R=self.null_points_2nd_order[0], Z=self.null_points_2nd_order[1]
+ )
+ self.Gbz_2nd_order = eq.tokamak.createBzGreensVec(
+ R=self.null_points_2nd_order[0], Z=self.null_points_2nd_order[1]
+ )
+ # for second order null
+ self.Gdbrdr_2nd_order = eq.tokamak.createdBrdrGreensVec(
+ R=self.null_points_2nd_order[0], Z=self.null_points_2nd_order[1]
+ )
+ self.Gdbzdz_2nd_order = eq.tokamak.createdBzdzGreensVec(
+ R=self.null_points_2nd_order[0], Z=self.null_points_2nd_order[1]
+ )
+ self.Gdbrdz_2nd_order = eq.tokamak.createdBrdzGreensVec(
+ R=self.null_points_2nd_order[0], Z=self.null_points_2nd_order[1]
+ )
+ self.Gdbzdr_2nd_order = eq.tokamak.createdBzdrGreensVec(
+ R=self.null_points_2nd_order[0], Z=self.null_points_2nd_order[1]
+ )
+
+ # ------------------------------------------------------------
+ # Flux value constraint Green's functions
+ # ------------------------------------------------------------
if self.psi_vals is not None:
# detect if psi constraints are defined on full grid
@@ -655,6 +700,89 @@ def build_plasma_vals(self, trial_plasma_psi):
/ self.null_points[0]
)
+ # ------------------------------------------------------------
+ # Magnetic field evaluation (and derivatives) at null points
+ #
+ # We evaluate the poloidal magnetic field components (B_R, B_Z)
+ # and their spatial derivatives at the identified null points
+ # using the Grad–Shafranov relations:
+ #
+ # B_R(R,Z) = - (1/R) ∂ψ/∂Z
+ # B_Z(R,Z) = (1/R) ∂ψ/∂R
+ #
+ # where ψ(R,Z) is the poloidal flux function.
+ #
+ # The callable `psi_func` provides ψ and its derivatives:
+ # dx = order of derivative with respect to R
+ # dy = order of derivative with respect to Z
+ #
+ # e.g.
+ # dx=1, dy=0 → ∂ψ/∂R
+ # dx=0, dy=1 → ∂ψ/∂Z
+ # dx=1, dy=1 → ∂²ψ/(∂R∂Z)
+ #
+ # These are then used to construct:
+ # - Magnetic field components (B_R, B_Z)
+ # - First spatial derivatives of the field:
+ # ∂B_R/∂R, ∂B_R/∂Z, ∂B_Z/∂R, ∂B_Z/∂Z
+ #
+ # All quantities are evaluated at the null point locations
+ # stored in `self.null_points_2nd_order = (R, Z)`.
+ # ------------------------------------------------------------
+ if self.null_points_2nd_order is not None:
+ dpsidr = psi_func(
+ self.null_points_2nd_order[0],
+ self.null_points_2nd_order[1],
+ dx=1,
+ dy=0,
+ grid=False,
+ )
+ dpsidz = psi_func(
+ self.null_points_2nd_order[0],
+ self.null_points_2nd_order[1],
+ dx=0,
+ dy=1,
+ grid=False,
+ )
+ d2psidrdz = psi_func(
+ self.null_points_2nd_order[0],
+ self.null_points_2nd_order[1],
+ dx=1,
+ dy=1,
+ grid=False,
+ )
+ d2psidr2 = psi_func(
+ self.null_points_2nd_order[0],
+ self.null_points_2nd_order[1],
+ dx=2,
+ dy=0,
+ grid=False,
+ )
+ d2psidz2 = psi_func(
+ self.null_points_2nd_order[0],
+ self.null_points_2nd_order[1],
+ dx=0,
+ dy=2,
+ grid=False,
+ )
+
+ # Br(R,Z)
+ self.Brp_2nd_order = -dpsidz / self.null_points_2nd_order[0]
+ # Bz(R,Z)
+ self.Bzp_2nd_order = dpsidr / self.null_points_2nd_order[0]
+ # dBrdr(R,Z)
+ self.dBrdrp_2nd_order = (
+ (dpsidz / self.null_points_2nd_order[0]) - d2psidrdz
+ ) / self.null_points_2nd_order[0]
+ # dBzdz(R,Z)
+ self.dBzdzp_2nd_order = d2psidrdz / self.null_points_2nd_order[0]
+ # dBrdz(R,Z)
+ self.dBrdzp_2nd_order = -d2psidz2 / self.null_points_2nd_order[0]
+ # dBzdr(R,Z)
+ self.dBzdrp_2nd_order = (
+ -(dpsidr / self.null_points_2nd_order[0]) + d2psidr2
+ ) / self.null_points_2nd_order[0]
+
# ------------------------------------------------------------
# Isoflux contour constraint evaluation
#
@@ -1022,6 +1150,16 @@ def build_lsq(self, full_currents_vec):
self.psiv_dim = len(b)
loss = loss + l
+ # second order null point constraints
+ if self.null_points_2nd_order is not None:
+ A_np_2nd_order, b_np_2nd_order, l = self.build_null_points_2nd_order_lsq(
+ full_currents_vec
+ )
+ A = np.concatenate((A, A_np_2nd_order), axis=0)
+ b = np.concatenate((b, b_np_2nd_order), axis=0)
+ self.nullp_2nd_order_dim = len(b)
+ loss = loss + l
+
# assemble the full system
self.A = np.copy(A)
self.b = np.copy(b)
@@ -1334,6 +1472,7 @@ def optimize_currents_grad(
trial_plasma_psi,
isoflux_weight=1.0,
null_points_weight=1.0,
+ null_points_2nd_order_weight=1.0,
psi_vals_weight=1.0,
):
"""
@@ -1398,6 +1537,11 @@ def optimize_currents_grad(
if self.null_points is not None:
b_weighted[idx : idx + self.nullp_dim] *= null_points_weight
idx += self.nullp_dim
+ if self.null_points_2nd_order is not None:
+ b_weighted[
+ idx : idx + self.nullp_2nd_order_dim
+ ] *= null_points_2nd_order_weight
+ idx += self.nullp_dim_2nd_order
if self.psi_vals is not None:
b_weighted[idx : idx + self.psiv_dim] *= psi_vals_weight
idx += self.psiv_dim
@@ -1624,6 +1768,128 @@ def build_plasma_isoflux_lsq(self, full_currents_vec, trial_plasma_psi):
self.b_plasma = np.concatenate(b, axis=0)
self.loss_plasma = np.linalg.norm(loss)
+ def build_null_points_2nd_order_lsq(self, full_currents_vec):
+ """
+ Construct the linear least-squares system enforcing second-order null point constraints.
+
+ This method assembles a linear system of the form:
+
+ A ΔI = b
+
+ where ΔI represents corrections to the *control coil currents*, such that
+ the magnetic field and its first spatial derivatives vanish (or match
+ desired values) at the specified second-order null points.
+
+ The constraints enforced at each null point are:
+
+ - B_R = 0
+ - B_Z = 0
+ - ∂B_R/∂R = 0
+ - ∂B_Z/∂Z = 0
+ - ∂B_R/∂Z = 0
+ - ∂B_Z/∂R = 0
+
+ Each constraint is written as a linear combination of coil currents using
+ precomputed Green's function matrices (G* terms), with an additional
+ contribution from the plasma.
+
+ The right-hand side vector `b` represents the *negative residual* of the
+ total field (coil + plasma) at the null points, so that solving the system
+ drives the residuals toward zero.
+
+ Parameters
+ ----------
+ full_currents_vec : np.ndarray
+ Array of shape (n_coils,) containing the full set of coil currents.
+ This includes both controlled and fixed coils, e.g. as returned by
+ `eq.tokamak.getCurrentsVec()`.
+
+ Returns
+ -------
+ A : np.ndarray
+ Least-squares matrix of shape (n_constraints, n_control_coils),
+ mapping control coil current corrections to changes in the magnetic
+ field quantities at the null points.
+
+ b : np.ndarray
+ Right-hand side vector of shape (n_constraints,), representing the
+ negative residuals of the current field configuration (coil + plasma).
+
+ loss : list of float
+ List of L2 norms of each individual constraint residual before solving.
+ The entries correspond to:
+ [||B_R||, ||B_Z||,
+ ||∂B_R/∂R||, ||∂B_Z/∂Z||,
+ ||∂B_R/∂Z||, ||∂B_Z/∂R||]
+
+ Notes
+ -----
+ - The matrices `G*` encode the linear response of the magnetic field (and
+ its derivatives) to coil currents.
+ - Contributions are split into:
+ * coil contribution: G * I
+ * plasma contribution: precomputed field values at null points
+ - Only coils selected by `self.control_mask` are treated as control variables.
+ - The system is typically overdetermined and intended to be solved in a
+ least-squares sense.
+ """
+
+ # Br field constraint
+ A_r = self.Gbr_2nd_order[self.control_mask].T
+ b_r = np.sum(
+ self.Gbr_2nd_order * full_currents_vec[:, np.newaxis], axis=0
+ ) # coils contribution
+ b_r += self.Brp_2nd_order # plasma contribution
+ loss = [np.linalg.norm(b_r)]
+
+ # Bz field constraint
+ A_z = self.Gbz_2nd_order[self.control_mask].T
+ b_z = np.sum(
+ self.Gbz_2nd_order * full_currents_vec[:, np.newaxis], axis=0
+ ) # coils contribution
+ b_z += self.Bzp_2nd_order # plasma contribution
+ loss.append(np.linalg.norm(b_z))
+
+ # dBrdr field constraint
+ A_r_deriv = self.Gdbrdr_2nd_order[self.control_mask].T
+ b_r_deriv = np.sum(
+ self.Gdbrdr_2nd_order * full_currents_vec[:, np.newaxis], axis=0
+ ) # coils contribution
+ b_r_deriv += self.dBrdrp_2nd_order # plasma contribution
+ loss.append(np.linalg.norm(b_r_deriv))
+
+ # dBzdz field constraint
+ A_z_deriv = self.Gdbzdz_2nd_order[self.control_mask].T
+ b_z_deriv = np.sum(
+ self.Gdbzdz_2nd_order * full_currents_vec[:, np.newaxis], axis=0
+ ) # coils contribution
+ b_z_deriv += self.dBzdzp_2nd_order # plasma contribution
+ loss.append(np.linalg.norm(b_z_deriv))
+
+ # dBrdz field constraint
+ A_r_deriv_cross = self.Gdbrdz_2nd_order[self.control_mask].T
+ b_r_deriv_cross = np.sum(
+ self.Gdbrdz_2nd_order * full_currents_vec[:, np.newaxis], axis=0
+ ) # coils contribution
+ b_r_deriv_cross += self.dBrdzp_2nd_order # plasma contribution
+ loss.append(np.linalg.norm(b_r_deriv_cross))
+
+ # dBzdr field constraint
+ A_z_deriv_cross = self.Gdbzdr_2nd_order[self.control_mask].T
+ b_z_deriv_cross = np.sum(
+ self.Gdbzdr_2nd_order * full_currents_vec[:, np.newaxis], axis=0
+ ) # coils contribution
+ b_z_deriv_cross += self.dBzdrp_2nd_order # plasma contribution
+ loss.append(np.linalg.norm(b_z_deriv_cross))
+
+ A = np.concatenate(
+ (A_r, A_z, A_r_deriv, A_z_deriv, A_r_deriv_cross, A_z_deriv_cross), axis=0
+ )
+ b = -np.concatenate(
+ (b_r, b_z, b_r_deriv, b_z_deriv, b_r_deriv_cross, b_z_deriv_cross), axis=0
+ )
+ return A, b, loss
+
def optimize_plasma_psi(self, full_currents_vec, trial_plasma_psi, l2_reg):
"""
Solve the regularised least-squares problem for plasma parameters.
diff --git a/requirements-freegs4e.txt b/requirements-freegs4e.txt
index db1fd7c1..275a1509 100644
--- a/requirements-freegs4e.txt
+++ b/requirements-freegs4e.txt
@@ -1 +1 @@
-freegs4e~=0.12
+freegs4e~=0.13