From 7bd00fad3c0ed68f9843c0fd7ce93babb9e99611 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 12 Feb 2026 17:34:06 +0000 Subject: [PATCH 01/13] enabling 2nd order null constraints --- freegsnke/inverse.py | 120 ++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 119 insertions(+), 1 deletion(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index 88c9b183..910974ff 100644 --- a/freegsnke/inverse.py +++ b/freegsnke/inverse.py @@ -33,6 +33,7 @@ def __init__( self, isoflux_set=None, null_points=None, + null_points_2nd_order=None, psi_vals=None, curr_vals=None, coil_current_limits=None, @@ -48,6 +49,10 @@ def __init__( null_points : list or np.array, optional structure [Rcoords, Zcoords], with Rcoords and Zcoords being 1D lists Sets the coordinates of the desired null points, including both Xpoints and Opoints + null_points_2nd_order : list or np.array, optional + structure [Rcoords, Zcoords], with Rcoords and Zcoords being 1D lists + Sets the coordinates of the desired 2nd order null points (typically used for making snowflake divertors). + Do not repeat these (R,Z) coords in the 'null_points' input, keep seperate. psi_vals : list or np.array, optional structure [Rcoords, Zcoords, psi_values] with Rcoords, Zcoords and psi_values having the same shape @@ -75,6 +80,10 @@ def __init__( if self.null_points is not None: self.null_points = np.array(self.null_points) + 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) + self.psi_vals = psi_vals if self.psi_vals is not None: self.full_grid = False @@ -208,6 +217,20 @@ def build_greens(self, eq): R=self.null_points[0], Z=self.null_points[1] ) + if self.null_points_2nd_order is not None: + 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] + ) + 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] + ) + if self.psi_vals is not None: if np.all(self.psi_vals[0] == eq.R.reshape(-1)) and np.all( self.psi_vals[1] == eq.Z.reshape(-1) @@ -221,7 +244,7 @@ def build_greens(self, eq): def build_plasma_vals(self, trial_plasma_psi): """Builds and stores all the values relative to the plasma, - based on the provided plasma_psi + based on the provided plasma_psi. Parameters ---------- @@ -233,6 +256,7 @@ def build_plasma_vals(self, trial_plasma_psi): self.eqR[:, 0], self.eqZ[0, :], trial_plasma_psi ) + # calculate values of constraint (for plasma flux) at 1st order null point locations if self.null_points is not None: self.brp = ( -psi_func(self.null_points[0], self.null_points[1], dy=1, grid=False) @@ -243,6 +267,40 @@ def build_plasma_vals(self, trial_plasma_psi): / self.null_points[0] ) + # calculate values of constraint (for plasma flux) at 2nd order null point locations + 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, + grid=False, + ) + dpsidz = psi_func( + self.null_points_2nd_order[0], + self.null_points_2nd_order[1], + 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, + ) + + # 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] + + # calculate flux values if self.isoflux_set is not None: self.d_psi_plasma_vals_iso = [] for i, isoflux in enumerate(self.isoflux_set): @@ -306,6 +364,52 @@ def build_null_points_lsq(self, full_currents_vec): b = -np.concatenate((b_r, b_z), axis=0) return A, b, loss + def build_null_points_2nd_order_lsq(self, full_currents_vec): + """Builds for the ordinary least sq problem associated to the 2nd order null points. + + Parameters + ---------- + full_currents_vec : np.array + Full vector of all coil current values. For example as returned by eq.tokamak.getCurrentsVec() + + """ + + # 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)) + + A = np.concatenate((A_r, A_z, A_r_deriv, A_z_deriv), axis=0) + b = -np.concatenate((b_r, b_z, b_r_deriv, b_z_deriv), axis=0) + return A, b, loss + def build_psi_vals_lsq(self, full_currents_vec): """Builds for the ordinary least sq problem associated to the psi values @@ -368,6 +472,14 @@ def build_lsq(self, full_currents_vec): b = np.concatenate((b, b_np), axis=0) self.nullp_dim = len(b) loss = loss + l + 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 if self.psi_vals is not None: A_pv, b_pv, l = self.build_psi_vals_lsq(full_currents_vec) A = np.concatenate((A, A_pv), axis=0) @@ -518,6 +630,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, current_weight=1.0, ): @@ -548,6 +661,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 From e714bd4948abc1414eaa8e3d077cba98a2b3afc2 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 12 Feb 2026 17:34:45 +0000 Subject: [PATCH 02/13] adding example showing how to use new constraint to make snowflake divertor --- ...ample01 - static_inverse_solve_MASTU.ipynb | 299 +++++++++++++++++- 1 file changed, 286 insertions(+), 13 deletions(-) diff --git a/examples/example01 - static_inverse_solve_MASTU.ipynb b/examples/example01 - static_inverse_solve_MASTU.ipynb index 0ed01ebc..b44855d5 100644 --- a/examples/example01 - static_inverse_solve_MASTU.ipynb +++ b/examples/example01 - static_inverse_solve_MASTU.ipynb @@ -10,11 +10,13 @@ "\n", "This example notebook shows how to use FreeGSNKE to solve **static inverse** free-boundary Grad-Shafranov (GS) problems. \n", "\n", - "In the **inverse** solve mode we seek to estimate the active poloidal field coil currents using user-defined constraints (e.g. on isoflux, x-point, o-point, and psi values) and plasma current density profiles for a desired equilibrium shape. \n", + "In the **inverse** solve mode we seek to estimate the active poloidal field coil currents using user-defined constraints (e.g. on isofluxes, magnetic null points, and psi values) and plasma current density profiles for a desired equilibrium shape.\n", "\n", "Note that during this solve, currents are **not** found in any specified passive structures. \n", "\n", - "Below, we illustrate how to use the solver for both diverted and limited plasma configurations in a **MAST-U-like tokamak** using stored pickle files containing the machine description. These machine description files partially come from the FreeGS repository and are not an exact replica of MAST-U." + "Below, we illustrate how to use the solver for both diverted and limited plasma configurations in a **MAST-U-like tokamak** using stored pickle files containing the machine description. These machine description files partially come from the FreeGS repository and are not an exact replica of MAST-U.\n", + "\n", + "We'll also try to construct more advanced configurations such as a \"snowflake\" divertor plasma. " ] }, { @@ -255,17 +257,81 @@ "source": [ "### Constraints\n", "\n", - "Recall that in the **inverse** solve mode we seek to **estimate the active poloidal field coil currents** using user-defined **constraints** (e.g. on isoflux, null points, and psi values) and plasma current density profiles for a desired equilibrium shape. \n", + "In **inverse solve** mode, we seek to **estimate the active poloidal field (PF) coil currents** using user-defined magnetic **constraints**, together with prescribed plasma current density profiles, in order to obtain a desired equilibrium shape.\n", + "\n", + "FreeGSNKE uses a `constrain` object to impose constraints on the magnetic configuration. The following constraint types are available:\n", + "\n", + "---\n", + "\n", + "#### First-order magnetic null points (`null_points`)\n", + "\n", + "These correspond to standard **X-points** or **O-points**. \n", + "\n", + "At each specified location $(R, Z)$, the solver enforces:\n", "\n", - "FreeGSNKE uses a `constrain` object, which accepts constraints on the location of:\n", - "- any null points, which are either X-points or O-points (`null_points`).\n", - "- any pairs of points that lie on the same flux surface (`isoflux_set`). The flux at these locations will be constrained to be the same. Multiple sets can be defined. \n", - "- any known flux points (`psi_vals`), i.e. one can specify $\\psi$ at any location $(R,Z)$, possibly an entire flux map (not used here). \n", - "- any known upper/lower limts on the coil currents being optimsied for (`coil_current_limits`).\n", + "$$\n", + "B_R(R, Z) = 0, \\qquad B_Z(R, Z) = 0.\n", + "$$\n", "\n", - "At least one constraint (preferably many more) is required to carry out an inverse solve.\n", + "---\n", + "\n", + "#### Second-order magnetic null points (`null_points_2nd_order`)\n", + "\n", + "These correspond to **second-order nulls of the poloidal magnetic field** (e.g. snowflake-type configurations).\n", + "\n", + "At each specified location $(R, Z)$, the solver enforces:\n", "\n", - "Here, we specify two X-point locations (we want a double null plasma), a number of isoflux locations, and some coil current limits. The isofluxes here will define the core plasma shape and part of the divertor legs. " + "$$\n", + "B_R(R, Z) = 0, \\qquad B_Z(R, Z) = 0,\n", + "$$\n", + "\n", + "$$\n", + "\\partial_R B_R(R, Z) = 0, \\qquad \\partial_Z B_Z(R, Z) = 0.\n", + "$$\n", + "\n", + "---\n", + "\n", + "#### Isoflux constraints (`isoflux_set`)\n", + "\n", + "These enforce that multiple points lie on the **same poloidal flux surface**.\n", + "\n", + "For a set of locations $(R_1, Z_1), (R_2, Z_2), \\dots, (R_N, Z_N)$, the solver enforces:\n", + "\n", + "$$\n", + "\\psi(R_1, Z_1) = \\psi(R_2, Z_2) = \\dots = \\psi(R_N, Z_N).\n", + "$$\n", + "\n", + "Multiple independent isoflux sets can be defined.\n", + "\n", + "---\n", + "\n", + "#### Prescribed poloidal flux values (`psi_vals`)\n", + "\n", + "These specify known values of the poloidal flux at selected locations:\n", + "\n", + "$$\n", + "\\psi(R, Z) = \\psi_{\\text{target}}.\n", + "$$\n", + "\n", + "This may also be used to prescribe a larger portion of the flux map (not used here).\n", + "\n", + "---\n", + "\n", + "#### Coil current limits (`coil_current_limits`)\n", + "\n", + "Upper and/or lower bounds may be imposed on the PF coil currents being optimised.\n", + "\n", + "---\n", + "\n", + "At least one constraint (and typically many more) is required to perform an inverse solve.\n", + "\n", + "In the example considered here, we specify:\n", + "\n", + "- Two X-point locations (corresponding to a **double-null plasma**),\n", + "- A number of isoflux constraints,\n", + "- Coil current limits.\n", + "\n", + "The isoflux constraints define the core plasma shape and part of the divertor legs.\n" ] }, { @@ -302,7 +368,6 @@ " [-5e3, -9e3, -9e3, -7e3, -7e3, -5e3, -4e3, -5e3, -10e3, -10e3, None]\n", "]\n", "\n", - " \n", "# instantiate the freegsnke constrain object\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", @@ -589,7 +654,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Second inverse solve: limited plasma example\n", + "### Second inverse solve: limited plasma\n", "\n", "Below we carry out an inverse solve seeking coil current values for a limited plasma configuration (rather than a diverted one). \n", "\n", @@ -604,7 +669,6 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "# a new eq object resets the coil currents and equilibrium\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak, # provide tokamak object\n", @@ -689,6 +753,215 @@ "plt.tight_layout()" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Third inverse solve: snowflake divertor\n", + "\n", + "Now we'll try to use the more complex constraints to form second order nulls and a snowflake divertor geometry. \n", + "\n", + "For this, we'll reset the equilibrium and profile objects (cranking up the resolution slightly)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 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.35 # outboard midplane radius\n", + "Rin = 0.32 # inboard midplane radius\n", + "\n", + "# locations of X-points (this time they will be second order nulls!)\n", + "Rx = 0.62\n", + "Zx = 1.2\n", + "\n", + "# second order null constraints to form the snowflake\n", + "null_points_2nd_order = [[Rx, Rx], [Zx, -Zx]]\n", + "\n", + "# isoflux constraints (no divertor constraints this time)\n", + "isoflux_set = np.array([[[Rx, Rx, Rin, Rout], [Zx, -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", + "# 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", + "# 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, # print output\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=True)\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=True)\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", + " [None] + coil_current_limits[0], # upper limits\n", + " [None] + coil_current_limits[1], # lower limits\n", + "):\n", + " # skip solenoid (no limits)\n", + " if ul is None or ll is None:\n", + " continue\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", + " if (coil_current >= ll) & (coil_current <= ul):\n", + " within_limit = True\n", + " else:\n", + " within_limit = False\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", + "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", "execution_count": null, From 95f71fe393f2e13952edd8ceac107d6a17e943ee Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 12 Feb 2026 17:36:34 +0000 Subject: [PATCH 03/13] slightly modified pickles --- .../data/simple_diverted_currents_PaxisIp.pk | Bin 5459 -> 5460 bytes .../data/simple_limited_currents_PaxisIp.pk | Bin 5459 -> 5460 bytes 2 files changed, 0 insertions(+), 0 deletions(-) diff --git a/examples/data/simple_diverted_currents_PaxisIp.pk b/examples/data/simple_diverted_currents_PaxisIp.pk index eae8a38d88303cd20604f0a18289b673ae7c5f49..08a31ddf02871b14f12ccd14449695a50c67cf65 100644 GIT binary patch delta 303 zcmcbtbwx|5fn}rH9GIa7qSC2A4C(n&9h&HnO|n;zm&MsXG5A-cj2K7dM8A%gkc_xb?;|xOf3n zJnd^rYg)nL12FSVq2g6)E-Vw&^6TLm%%Bp_WE%qCtuKaaFo%lop8s`=*+;ng028RV scy-T~LpN5#)Ca&UFk)bTTk&KaTpVV8ywxkm=Nw<`AP(Bx#duN}0Pmf5L;wH) delta 302 zcmcbjby-WXfn}(flQcmY)W z;;i)2{u31kVCI`b#ivgZ{Lsj|rXH@r3@Wi_cfpk0yjVSoE0?c_i^I(S^zW*v%4Kj3=1_6bIq$PM+dE52`ZVVNF9>{8&&ov7!UH}!( z(R$>1XnNuSnE9qq@s?HGU-VARu7_(dgGx;NyJOG($9-@O=1}pNXpx*PD_ddW0VYs! r-pcJPFPv7w#bFlY?!LY9z1|YIIL!P)*UMYEV^a4+9JINc@uV;S!*hBA From 0effbd7d2cb3c487e9fdd6db8f20008425d11534 Mon Sep 17 00:00:00 2001 From: kpentland Date: Fri, 13 Feb 2026 10:54:02 +0000 Subject: [PATCH 04/13] tweaking notebook example for snowflake --- ...ample01 - static_inverse_solve_MASTU.ipynb | 114 +++++++++++++----- 1 file changed, 83 insertions(+), 31 deletions(-) diff --git a/examples/example01 - static_inverse_solve_MASTU.ipynb b/examples/example01 - static_inverse_solve_MASTU.ipynb index b44855d5..1344a911 100644 --- a/examples/example01 - static_inverse_solve_MASTU.ipynb +++ b/examples/example01 - static_inverse_solve_MASTU.ipynb @@ -362,17 +362,17 @@ "coil_current_limits = [\n", "# upper limits...\n", "# PX, D1, D2, D3, Dp, D5, D6, D7, P4, P5, P6\n", - " [5e3, 9e3, 9e3, 7e3, 7e3, 5e3, 4e3, 5e3, 0, 0, None],\n", + " [6e3, 9e3, 9e3, 7e3, 7e3, 5e3, 4e3, 5e3, 0, 0, None],\n", "# lower limits...\n", "# PX, D1, D2, D3, Dp, D5, D6, D7, P4, P5, P6\n", - " [-5e3, -9e3, -9e3, -7e3, -7e3, -5e3, -4e3, -5e3, -10e3, -10e3, None]\n", + " [-6e3, -9e3, -9e3, -7e3, -7e3, -5e3, -4e3, -5e3, -12e3, -12e3, None]\n", "]\n", "\n", "# instantiate the freegsnke constrain object\n", "constrain = Inverse_optimizer(\n", " null_points=null_points,\n", " isoflux_set=isoflux_set,\n", - " coil_current_limits=coil_current_limits\n", + " coil_current_limits=coil_current_limits,\n", ")\n", "\n", "# if you find coil limits are being violated or an undesireable solution is being produced,\n", @@ -589,19 +589,19 @@ " [None] + coil_current_limits[0], # upper limits\n", " [None] + coil_current_limits[1], # lower limits\n", "):\n", - " # skip solenoid (no limits)\n", - " if ul is None or ll is None:\n", - " continue\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", - " if (coil_current >= ll) & (coil_current <= ul):\n", - " within_limit = True\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 = False\n", + " within_limit = None\n", "\n", " # print\n", " print(f\"{coil_name}: {coil_current:.2f} --> [{ll},{ul}] [A] --> limits met = {within_limit}\")\n", @@ -616,7 +616,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can also visualise the currents with respect to their limits (normalising between -1 and 1). Note the Solenoid is missing as it has no limits. " + "We can also visualise the currents with respect to their limits (normalising between -1 and 1). Note the Solenoid and P6 are missing as they have no limits. " ] }, { @@ -628,6 +628,8 @@ "# 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", @@ -770,6 +772,14 @@ "metadata": {}, "outputs": [], "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", @@ -793,27 +803,36 @@ "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# first we specify some alternative constraints\n", - "Rout = 1.35 # outboard midplane radius\n", - "Rin = 0.32 # inboard midplane radius\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.62\n", + "Rx = 0.65\n", "Zx = 1.2\n", "\n", "# second order null constraints to form the snowflake\n", - "null_points_2nd_order = [[Rx, Rx], [Zx, -Zx]]\n", + "null_points_2nd_order = [[Rx], [Zx]]\n", "\n", "# isoflux constraints (no divertor constraints this time)\n", - "isoflux_set = np.array([[[Rx, Rx, Rin, Rout], [Zx, -Zx, 0.,0.]]])\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", + "# 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", @@ -821,6 +840,9 @@ " 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", @@ -867,7 +889,7 @@ "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", + "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", @@ -877,7 +899,7 @@ "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", + "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()" @@ -905,22 +927,22 @@ "for coil_name, coil_current, ul, ll in zip(\n", " eq.tokamak.coils_dict,\n", " eq.tokamak.getCurrentsVec()[:12],\n", - " [None] + coil_current_limits[0], # upper limits\n", - " [None] + coil_current_limits[1], # lower limits\n", + " coil_current_limits[0] + [None], # upper limits\n", + " coil_current_limits[1] + [None], # lower limits\n", "):\n", - " # skip solenoid (no limits)\n", - " if ul is None or ll is None:\n", - " continue\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", - " if (coil_current >= ll) & (coil_current <= ul):\n", - " within_limit = True\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 = False\n", + " within_limit = None\n", "\n", " # print\n", " print(f\"{coil_name}: {coil_current:.2f} --> [{ll},{ul}] [A] --> limits met = {within_limit}\")\n", @@ -940,6 +962,8 @@ "# 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", @@ -969,6 +993,34 @@ "outputs": [], "source": [] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "code", "execution_count": null, From 9f7091a32c1d4521048f64d87a5192a5a2113563 Mon Sep 17 00:00:00 2001 From: kpentland Date: Fri, 13 Feb 2026 10:54:37 +0000 Subject: [PATCH 05/13] adding missing cross derivative constraints --- freegsnke/inverse.py | 54 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 52 insertions(+), 2 deletions(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index 910974ff..5d06f394 100644 --- a/freegsnke/inverse.py +++ b/freegsnke/inverse.py @@ -218,18 +218,26 @@ def build_greens(self, eq): ) 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] + ) if self.psi_vals is not None: if np.all(self.psi_vals[0] == eq.R.reshape(-1)) and np.all( @@ -273,11 +281,13 @@ def build_plasma_vals(self, trial_plasma_psi): 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, ) @@ -288,6 +298,20 @@ def build_plasma_vals(self, trial_plasma_psi): 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] @@ -299,6 +323,12 @@ def build_plasma_vals(self, trial_plasma_psi): ) / 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] # calculate flux values if self.isoflux_set is not None: @@ -406,8 +436,28 @@ def build_null_points_2nd_order_lsq(self, full_currents_vec): b_z_deriv += self.dBzdzp_2nd_order # plasma contribution loss.append(np.linalg.norm(b_z_deriv)) - A = np.concatenate((A_r, A_z, A_r_deriv, A_z_deriv), axis=0) - b = -np.concatenate((b_r, b_z, b_r_deriv, b_z_deriv), axis=0) + # 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 build_psi_vals_lsq(self, full_currents_vec): From dec7a9012f2b7f3a32aa0048360e3134a52ac9fb Mon Sep 17 00:00:00 2001 From: kpentland Date: Tue, 24 Mar 2026 11:22:15 +0000 Subject: [PATCH 06/13] formatting fix --- freegsnke/inverse.py | 16 +++++++--------- 1 file changed, 7 insertions(+), 9 deletions(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index 47cecda1..f7cd7f31 100644 --- a/freegsnke/inverse.py +++ b/freegsnke/inverse.py @@ -211,7 +211,7 @@ def __init__( 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 @@ -588,7 +588,6 @@ def build_greens(self, eq): R=self.null_points[0], Z=self.null_points[1] ) - # ------------------------------------------------------------ # Magnetic null point Green's functions (for second-order nulls) # ------------------------------------------------------------ @@ -614,7 +613,7 @@ def build_greens(self, eq): 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 # ------------------------------------------------------------ @@ -789,7 +788,7 @@ def build_plasma_vals(self, trial_plasma_psi): self.dBzdrp_2nd_order = ( -(dpsidr / self.null_points_2nd_order[0]) + d2psidr2 ) / self.null_points_2nd_order[0] - + # ------------------------------------------------------------ # Isoflux contour constraint evaluation # @@ -997,7 +996,6 @@ def build_null_points_lsq(self, full_currents_vec): return A, b, loss - def build_psi_vals_lsq(self, full_currents_vec): """ Construct a least-squares system enforcing direct flux value constraints. @@ -1269,7 +1267,7 @@ def build_lsq(self, full_currents_vec): b = np.concatenate((b, b_pv), axis=0) * self.weight_psi self.psiv_dim = len(b) loss = loss + l - + # direct flux value 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( @@ -1279,7 +1277,7 @@ def build_lsq(self, full_currents_vec): 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) @@ -1887,7 +1885,7 @@ def build_plasma_isoflux_lsq(self, full_currents_vec, trial_plasma_psi): self.A_plasma = np.concatenate(A, axis=0) 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. @@ -2009,7 +2007,7 @@ def build_null_points_2nd_order_lsq(self, full_currents_vec): (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. From b896f7f30661578f15dd51899f27591ad65b18a3 Mon Sep 17 00:00:00 2001 From: kpentland Date: Tue, 24 Mar 2026 11:39:37 +0000 Subject: [PATCH 07/13] adding in snowflake divertor example to 01b --- ...e01b - advanced_static_inverse_solve.ipynb | 275 +++++++++++++++++- 1 file changed, 272 insertions(+), 3 deletions(-) diff --git a/examples/example01b - advanced_static_inverse_solve.ipynb b/examples/example01b - advanced_static_inverse_solve.ipynb index 95f6b332..3d1f72f0 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,281 @@ "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", + "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", From 78bea777a9276351e6b66cd8c462d0578bee7ee1 Mon Sep 17 00:00:00 2001 From: kpentland Date: Mon, 27 Apr 2026 09:11:23 +0100 Subject: [PATCH 08/13] updating FreeGS4E release required --- requirements-freegs4e.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From c1f4fbc13e876664560a58b1475acda2261a127f Mon Sep 17 00:00:00 2001 From: kpentland Date: Mon, 27 Apr 2026 09:13:54 +0100 Subject: [PATCH 09/13] updating example01b notebook wrt review comments --- examples/example01b - advanced_static_inverse_solve.ipynb | 2 ++ 1 file changed, 2 insertions(+) diff --git a/examples/example01b - advanced_static_inverse_solve.ipynb b/examples/example01b - advanced_static_inverse_solve.ipynb index 3d1f72f0..1a75390c 100644 --- a/examples/example01b - advanced_static_inverse_solve.ipynb +++ b/examples/example01b - advanced_static_inverse_solve.ipynb @@ -464,6 +464,8 @@ "\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)." ] }, From 8c20e6c607bd21f9eeeb3af4a4b94b1edecf94a9 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 14 May 2026 10:24:05 +0100 Subject: [PATCH 10/13] fixing broken docstrings --- freegsnke/inverse.py | 72 ++++++++++++++++++++------------------------ 1 file changed, 33 insertions(+), 39 deletions(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index f7cd7f31..881a3e34 100644 --- a/freegsnke/inverse.py +++ b/freegsnke/inverse.py @@ -74,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. @@ -90,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. @@ -101,37 +112,16 @@ def __init__( Used to enforce ψ(R,Z) = ψ_target at specified locations. - isoflux_set : list or np.array, optional - list of isoflux objects, each with structure - [Rcoords, Zcoords] - with Rcoords and Zcoords being 1D lists of the coords of all points that are requested to have the same flux value - null_points : list or np.array, optional - structure [Rcoords, Zcoords], with Rcoords and Zcoords being 1D lists - Sets the coordinates of the desired null points, including both Xpoints and Opoints - null_points_2nd_order : list or np.array, optional - structure [Rcoords, Zcoords], with Rcoords and Zcoords being 1D lists - Sets the coordinates of the desired 2nd order null points (typically used for making snowflake divertors). - Do not repeat these (R,Z) coords in the 'null_points' input, keep seperate. - psi_vals : list or np.array, optional - structure [Rcoords, Zcoords, psi_values] - with Rcoords, Zcoords and psi_values having the same shape - Sets the desired values of psi for a set of coordinates, possibly an entire map - curr_vals : list, optional - structure [[coil indexes in the array of coils available for control], [coil current values]] coil_current_limits : list, optional Hard inequality bounds on coil currents. 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. @@ -142,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 # ------------------------------------------------------------ From 394ffcbf2f93c762fcc268c54ef54fce11f7c772 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 14 May 2026 10:43:11 +0100 Subject: [PATCH 11/13] black --- freegsnke/inverse.py | 233 ++++++++++++++++++++++--------------------- 1 file changed, 117 insertions(+), 116 deletions(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index 881a3e34..b1a280d2 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 @@ -990,118 +990,118 @@ def build_null_points_lsq(self, full_currents_vec): return A, b, loss - def build_psi_vals_lsq(self, full_currents_vec): - """ - Construct a least-squares system enforcing direct flux value constraints. - - This method builds the optimisation system: - - A I_control ≈ b - - where: - - ψ_model(R_i, Z_i) = ψ_target(R_i, Z_i) - - is enforced by matching magnetic flux values at specified locations. - - The optimisation residual is defined as: - - b = ψ_tokamak(I) + ψ_plasma - - ψ_target - - ⟨b⟩ - - Mean flux removal is applied to remove arbitrary vertical flux offsets, - since the Grad–Shafranov equation is invariant under constant flux shifts. - - Parameters - ---------- - full_currents_vec : ndarray - Full coil current vector. - - Example: - eq.tokamak.getCurrentsVec() - - Returns - ------- - A : ndarray - Jacobian matrix mapping coil current perturbations → flux changes. - - b : ndarray - Flux mismatch residual vector. - - normalised_loss : list of float - Normalised constraint violation magnitude. - - Notes - ----- - This constraint formulation is commonly used for: - - • Magnetic axis pinning - • Boundary flux matching - • Profile shape control - - Mathematical formulation - ------------------------ - Solve: - - min_I || G I + ψ_plasma − ψ_target ||² - """ - - # flux response wrt coil currents - 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 build_psi_vals_lsq(self, full_currents_vec): + # """ + # Construct a least-squares system enforcing direct flux value constraints. + + # This method builds the optimisation system: + + # A I_control ≈ b + + # where: + + # ψ_model(R_i, Z_i) = ψ_target(R_i, Z_i) + + # is enforced by matching magnetic flux values at specified locations. + + # The optimisation residual is defined as: + + # b = ψ_tokamak(I) + ψ_plasma + # - ψ_target + # - ⟨b⟩ + + # Mean flux removal is applied to remove arbitrary vertical flux offsets, + # since the Grad–Shafranov equation is invariant under constant flux shifts. + + # Parameters + # ---------- + # full_currents_vec : ndarray + # Full coil current vector. + + # Example: + # eq.tokamak.getCurrentsVec() + + # Returns + # ------- + # A : ndarray + # Jacobian matrix mapping coil current perturbations → flux changes. + + # b : ndarray + # Flux mismatch residual vector. + + # normalised_loss : list of float + # Normalised constraint violation magnitude. + + # Notes + # ----- + # This constraint formulation is commonly used for: + + # • Magnetic axis pinning + # • Boundary flux matching + # • Profile shape control + + # Mathematical formulation + # ------------------------ + # Solve: + + # min_I || G I + ψ_plasma − ψ_target ||² + # """ + + # # flux response wrt coil currents + # 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 build_psi_vals_lsq(self, full_currents_vec): """ @@ -1254,7 +1254,8 @@ def build_lsq(self, full_currents_vec): b = np.concatenate((b, b_np), axis=0) * self.weight_nulls self.nullp_dim = len(b) loss = loss + l - # second order null point constraints + + # direct flux value constraints if self.psi_vals is not None: A_pv, b_pv, l = self.build_psi_vals_lsq(full_currents_vec) A = np.concatenate((A, A_pv), axis=0) @@ -1262,7 +1263,7 @@ def build_lsq(self, full_currents_vec): self.psiv_dim = len(b) loss = loss + l - # direct flux value constraints + # 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 From e45cf0a830f00842316682fad52e94895cae1156 Mon Sep 17 00:00:00 2001 From: kpentland Date: Fri, 15 May 2026 08:55:05 +0100 Subject: [PATCH 12/13] Removing old commented out code from inverse.py --- freegsnke/inverse.py | 113 ------------------------------------------- 1 file changed, 113 deletions(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index b1a280d2..33a8e851 100644 --- a/freegsnke/inverse.py +++ b/freegsnke/inverse.py @@ -990,119 +990,6 @@ def build_null_points_lsq(self, full_currents_vec): return A, b, loss - # def build_psi_vals_lsq(self, full_currents_vec): - # """ - # Construct a least-squares system enforcing direct flux value constraints. - - # This method builds the optimisation system: - - # A I_control ≈ b - - # where: - - # ψ_model(R_i, Z_i) = ψ_target(R_i, Z_i) - - # is enforced by matching magnetic flux values at specified locations. - - # The optimisation residual is defined as: - - # b = ψ_tokamak(I) + ψ_plasma - # - ψ_target - # - ⟨b⟩ - - # Mean flux removal is applied to remove arbitrary vertical flux offsets, - # since the Grad–Shafranov equation is invariant under constant flux shifts. - - # Parameters - # ---------- - # full_currents_vec : ndarray - # Full coil current vector. - - # Example: - # eq.tokamak.getCurrentsVec() - - # Returns - # ------- - # A : ndarray - # Jacobian matrix mapping coil current perturbations → flux changes. - - # b : ndarray - # Flux mismatch residual vector. - - # normalised_loss : list of float - # Normalised constraint violation magnitude. - - # Notes - # ----- - # This constraint formulation is commonly used for: - - # • Magnetic axis pinning - # • Boundary flux matching - # • Profile shape control - - # Mathematical formulation - # ------------------------ - # Solve: - - # min_I || G I + ψ_plasma − ψ_target ||² - # """ - - # # flux response wrt coil currents - # 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 build_psi_vals_lsq(self, full_currents_vec): """ Construct a least-squares system enforcing direct flux value constraints. From 2e2b2fffedb96d67750103e2b4786f95cae42ebb Mon Sep 17 00:00:00 2001 From: kpentland Date: Fri, 15 May 2026 11:04:40 +0100 Subject: [PATCH 13/13] Minor commenting change to re-run CI-CD pipeline --- freegsnke/inverse.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/freegsnke/inverse.py b/freegsnke/inverse.py index 33a8e851..2866920c 100644 --- a/freegsnke/inverse.py +++ b/freegsnke/inverse.py @@ -28,7 +28,7 @@ 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__(