Merge branch 'main' into share-profile-geometry

Author: kpentland

.github/workflows/notebooks.yml

@@ -30,11 +30,11 @@ jobs:
         run: jupyter nbconvert --execute --to notebook --inplace "examples/example03 - extracting_equilibrium_quantities.ipynb"
       - name: Check using magnetic probes notebook
         run: jupyter nbconvert --execute --to notebook --inplace "examples/example04 - using_magnetic_probes.ipynb"
-      - name: Check evolutive forward solve notebook
+      - name: Check evolutive forward solve notebook (5a)
         run: jupyter nbconvert --execute --to notebook --inplace "examples/example05a - nonlinear_and_linear_evolution_with_GS.ipynb"
-      - name: Check evolutive forward solve notebook
+      - name: Check evolutive forward solve notebook (5b)
         run: jupyter nbconvert --execute --to notebook --inplace "examples/example05b - linear_evolution_without_GS.ipynb"
-      - name: Check evolutive forward solve notebook
+      - name: Check evolutive forward solve notebook (5c)
         run: jupyter nbconvert --execute --to notebook --inplace "examples/example05c - linear_evolution_with_relinearisation.ipynb"
       - name: Check static inverse solve Anamak notebook
         run: jupyter nbconvert --execute --to notebook --inplace "examples/example07a - Anamak"

README.md

@@ -29,7 +29,7 @@ These problems can be solved in a **user-specified tokamak geometry** that can i
 | ------ | ------ | ------ | ------ |
 | Active poloidal field coils | Can be assigned (voltage-driven) currents that influence plasma shape and position. | Locations, sizes (areas), wirings (series/anti-series), polarities (+1 or -1), resistivities (of coil materials), and number of windings. | Blue rectangles |
 | Passive conducting structures  | Can be assigned induced eddy currents that also impact plasma shape and position. In evolutive forward mode, these are solved self-consistently. | Locations, sizes, orientations (if available), and filaments (as passives can be refined if needed). | Dark grey parallelograms |
-| Wall and/or limiter contours  | Confines the plasma boundary (for computational purposes). | Locations. | Solid black line |
+| Wall and/or limiter contours  | Confines the plasma boundary (for computational purposes). | Locations; the limiter contour must lie strictly inside the equilibrium solution domain. | Solid black line |
 | Magnetic diagnostic probes  | Can measure the poloidal flux (fluxloops) or the magnetic field strength (pickup coils) at specified locations. | Locations (for both) and orientations (for pickup coils). | Orange diamonds (fluxloops) and brown dots/lines (pickup coils) |
 
 Static Grad-Shafranov problems are solved using **fourth-order accurate finite differences** and a **purpose-built Newton-Krylov method** for additional **stability and convergence** speed (over the Picard iterations used in FreeGS). An implicit Euler method and the same Newton-Krylov solver are used to tackle the evolutive problem.

examples/README.md

@@ -21,4 +21,4 @@ These example Jupyter notebooks are intended to be the **first port of call for
 | Example 10 | Learn how to calculate growth rates associated with vertically unstable modes. | Anyone |
 | Example 11 | Learn how to use the evolutive solver alongside a virtual plasma control system (FreeGSNKE Pulse Design Tool). | Anyone |
 
-If a new example has been created, please add a new line to the table explaining its purpose!
+If a new example has been created, please add a new line to the table explaining its purpose!

examples/example00 - build_tokamak_machine.ipynb

@@ -375,7 +375,7 @@
     "\n",
     "For diverted plasmas, the limiter has no effect on the equilibrium or its evolution. For limiter plasmas, the last closed flux surface will be tangent to the limiter. FreeGSNKE will also force the last closed flux surface to lie entirely inside the region allowed by the limiter. In other words, outside the limiter, the constraint $J=0$ on the plasma current density is enforced. \n",
     "\n",
-    "The format of the limiter coordinates for FreeGSNKE is a simple list of dictionaries with `R` and `Z` coordinates."
+    "The format of the limiter coordinates for FreeGSNKE is a simple list of dictionaries with `R` and `Z` coordinates. When an `Equilibrium` is later instantiated, the whole limiter contour must lie strictly inside the rectangular solution domain defined by `Rmin`, `Rmax`, `Zmin`, and `Zmax`. Active coils and passive structures may sit outside this grid, but the limiter itself may not."
    ]
   },
   {
@@ -699,6 +699,161 @@
     "ax1.set_xlabel(r'Major radius, $R$ [m]')\n",
     "ax1.set_ylabel(r'Height, $Z$ [m]')"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Updating an existing machine description\n",
+    "\n",
+    "If you already have a `tokamak` object, you can update its machine description directly from dictionaries or pickle paths.\n",
+    "\n",
+    "This is useful when you want to test small geometry changes without writing new pickle files or rebuilding every object from scratch.\n",
+    "\n",
+    "If the `tokamak` is already attached to an `Equilibrium` object (see the next example notebook), use `eq.update_machine_description(...)` instead so the equilibrium refreshes its cached Green's functions and limiter masks automatically. Note that after changing the machine geometry, any existing solver objects should be reinstantiated (again, see next notebooks for solver objects!).\n",
+    "\n",
+    "Let's start by moving the `Pz` coil radially outward!"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from copy import deepcopy\n",
+    "\n",
+    "# preserve an existing current so we can see it survive the machine update\n",
+    "tokamak_alt.set_coil_current(\"Pz\", 1.0e3)\n",
+    "\n",
+    "# extract the old R coords for Pz\n",
+    "old_pz_R = tokamak_alt.coils_dict[\"Pz\"][\"coords\"][0].copy()\n",
+    "\n",
+    "# work on a copy of the machine-description data already built above.\n",
+    "updated_active_coils = deepcopy(active_coils)\n",
+    "\n",
+    "# modify the Pz coords in this dictionary\n",
+    "updated_active_coils[\"Pz\"][\"\"][\"R\"] = [\n",
+    "    r + 0.6 for r in updated_active_coils[\"Pz\"][\"\"][\"R\"]\n",
+    "]\n",
+    "\n",
+    "# call this function (faster than re-building the whole machine as it only looks at changed coils)\n",
+    "tokamak_alt.set_machine_description(\n",
+    "    active_coils_data=updated_active_coils,\n",
+    "    passive_coils_data=passive_coils,\n",
+    "    limiter_data=limiter,\n",
+    "    wall_data=wall,\n",
+    "    magnetic_probe_data=magnetic_probes,\n",
+    ")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# check that the coil has been moved\n",
+    "new_pz_R = tokamak_alt.coils_dict[\"Pz\"][\"coords\"][0]\n",
+    "print(f\"Pz moved from R = {old_pz_R} to R = {new_pz_R}\")\n",
+    "\n",
+    "# check the coil current didn't change\n",
+    "print(f\"Preserved Pz current: {tokamak_alt['Pz'].current} A\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot the updated machine and notice Pz outside the vessel\n",
+    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
+    "plt.tight_layout()\n",
+    "\n",
+    "tokamak_alt.plot(axis=ax1, show=False)\n",
+    "ax1.plot(tokamak_alt.limiter.R, tokamak_alt.limiter.Z, color=\"k\", linewidth=1.2, linestyle=\"--\")\n",
+    "ax1.plot(tokamak_alt.wall.R, tokamak_alt.wall.Z, color=\"k\", linewidth=1.2, linestyle=\"-\")\n",
+    "\n",
+    "ax1.grid(alpha=0.5)\n",
+    "ax1.set_aspect(\"equal\")\n",
+    "ax1.set_xlabel(r\"Major radius, $R$ [m]\")\n",
+    "ax1.set_ylabel(r\"Height, $Z$ [m]\")\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Updating an active coil dictionary individually\n",
+    "\n",
+    "For optimization workflows that perturb one active coil (or circuit) at a time, we can use `tokamak.update_active_coil(coil_name, active_coil_data)`. \n",
+    "\n",
+    "This takes the data stored under that coil's label in the active-coils dictionary, updates only that active coil/circuit, and recalculates only the affected R/M matrix entries. On an equilibrium object, use `eq.update_active_coil(...)` so only the affected coil Greens function is refreshed (see next notebooks!).\n",
+    "\n",
+    "For example, you may want to do this in **machine geometry optimisation** studies, e.g. [Nunn et al. (2025)](https://pubs.aip.org/aip/pop/article/32/7/072507/3355562/Bayesian-optimization-of-poloidal-field-coil).\n",
+    "\n",
+    "Let's move the `P1` circuit vertically outward this time."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# old P1 coil data\n",
+    "old_p1_Z = tokamak_alt.coils_dict[\"P1\"][\"coords\"][1].copy()\n",
+    "\n",
+    "# extract only the P1 circuit data\n",
+    "updated_p1 = deepcopy(updated_active_coils[\"P1\"])\n",
+    "\n",
+    "# make changes\n",
+    "updated_p1[\"upper\"][\"Z\"] = [z + 0.4 for z in updated_p1[\"upper\"][\"Z\"]]\n",
+    "updated_p1[\"lower\"][\"Z\"] = [z - 0.4 for z in updated_p1[\"lower\"][\"Z\"]]\n",
+    "\n",
+    "# update the circuit in the tokamak\n",
+    "tokamak_alt.update_active_coil(\"P1\", updated_p1)\n",
+    "\n",
+    "# show the changes\n",
+    "new_p1_Z = tokamak_alt.coils_dict[\"P1\"][\"coords\"][1]\n",
+    "print(f\"P1 moved from Z = {old_p1_Z} to Z = {new_p1_Z}\")\n",
+    "print(f\"Changed coils recorded: {tokamak_alt._last_machine_update_changed_coils}\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot the updated machine\n",
+    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
+    "plt.tight_layout()\n",
+    "\n",
+    "tokamak_alt.plot(axis=ax1, show=False)\n",
+    "ax1.plot(tokamak_alt.limiter.R, tokamak_alt.limiter.Z, color=\"k\", linewidth=1.2, linestyle=\"--\")\n",
+    "ax1.plot(tokamak_alt.wall.R, tokamak_alt.wall.Z, color=\"k\", linewidth=1.2, linestyle=\"-\")\n",
+    "\n",
+    "ax1.grid(alpha=0.5)\n",
+    "ax1.set_aspect(\"equal\")\n",
+    "ax1.set_xlabel(r\"Major radius, $R$ [m]\")\n",
+    "ax1.set_ylabel(r\"Height, $Z$ [m]\")\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
   }
  ],
  "metadata": {
@@ -717,7 +872,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.10.16"
+   "version": "3.10.20"
   }
  },
  "nbformat": 4,

examples/example01a - static_inverse_solve_MASTU.ipynb

@@ -116,7 +116,7 @@
                 "\n",
                 "We are now ready to build a plasma equilibrium object for our tokamak. This is done using the `freegs4e.Equilibrium` class, which implicitly defines the rectangular domain of the solver as well as the grid resolution.\n",
                 "\n",
-                "`Equilibrium` has sensible defaults, but it is recommended to define the radial and vertical domain of the grid using the `Rmin`, `Rmax`, `Zmin` and `Zmax` parameters, as well as the grid resolution in the radial and vertical directions with the `nx` and `ny` parameters. The grid will be initialised using fourth-order finite differences. Note that the computational grid **must** encompass the domain enclosed by the limiter object (as this is where the plasma will be confined to). Though it does not need to encompass any/all of the active coils or passive structures. \n",
+                "`Equilibrium` has sensible defaults, but it is recommended to define the radial and vertical domain of the grid using the `Rmin`, `Rmax`, `Zmin` and `Zmax` parameters, as well as the grid resolution in the radial and vertical directions with the `nx` and `ny` parameters. The grid will be initialised using fourth-order finite differences. Note that the limiter contour **must** lie strictly inside the computational grid, because limiter-boundary flux interpolation is used when identifying limited plasmas. The grid does not need to encompass any/all of the active coils or passive structures. \n",
                 "\n",
                 "A tokamak object should be supplied to the `tokamak` parameter to assign the desired machine to the equilibrium.\n",
                 "\n",