Merge branch 'main' into coil_activation

Author: kpentland

README.md

@@ -186,8 +186,7 @@ Here are a list of FreeGSNKE papers that describe or use the code:
 - N. C. Amorisco et al, "FreeGSNKE: A Python-based dynamic free-boundary toroidal plasma equilibrium solver", Physics of Plasmas, **31**, 042517 (2024). DOI: [10.1063/5.0188467](https://doi.org/10.1063/5.0188467).
 - A. Agnello et al, "Emulation techniques for scenario and classical control design of tokamak plasmas", Physics of Plasmas, **31**, 043091 (2024). DOI: [10.1063/5.0187822](https://doi.org/10.1063/5.0187822).
 - K. Pentland et al, "Validation of the static forward Grad-Shafranov equilibrium solvers in FreeGSNKE and Fiesta using EFIT++ reconstructions from MAST-U", Physica Scripta, **100**, 025608 (2025). DOI: [10.1088/1402-4896/ada192](https://iopscience.iop.org/article/10.1088/1402-4896/ada192).
-- K. Pentland et al, "Multiple solutions to the static forward free-boundary Grad-Shafranov problem on MAST-U", arXiv (2025). DOI: [2503.05674](https://arxiv.org/abs/2503.05674).
- 
+- K. Pentland et al, "Multiple solutions to the static forward free-boundary Grad-Shafranov problem on MAST-U", Nuclear Fusion (2025). DOI: [10.1088/1741-4326/adf3cc](https://iopscience.iop.org/article/10.1088/1741-4326/adf3cc). 
 
 If you would like your FreeGSNKE-related paper to be added, please let us know!
 

examples/example01 - static_inverse_solve_MASTU.ipynb

@@ -151,14 +151,18 @@
     "\n",
     "We can now instatiate a profile object that contains the chosen parameterisation of the toroidal plasma current density $J_p$ (i.e. on right hand side of the GS equation). We can then set the paramters for the chosen current density profiles. \n",
     "\n",
-    "A number of commonly used profile parameterisations exist in FreeGSNKE, including:\n",
-    "- `ConstrainPaxisIp`\n",
-    "- `ConstrainBetapIp`\n",
-    "- `Fiesta_Topeol`\n",
-    "- `Lao85`\n",
-    "- `TensionSpline`\n",
-    "\n",
-    "In this notebook, we will make use of the `ConstrainPaxisIp` (and `ConstrainBetapIp`) profiles (see [Jeon (2015)](https://link.springer.com/article/10.3938/jkps.67.843)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue. \n",
+    "The following table indicates which parameterisations are currently available in FreeGSNKE and for which types of equilibrium simulation they have been tested to work:\n",
+    "\n",
+    "| Name | Static | Evolutive (linear) | Evolutive (nonlinear) |\n",
+    "|-------|--------|------------------|---------------------|\n",
+    "| `ConstrainPaxisIp` | ✅ | ✅ | ✅ |\n",
+    "| `ConstrainBetapIp` | ✅ | ✅ | ✅ |\n",
+    "| `FiestaTopeol` | ✅ | ✅ | ✅ |\n",
+    "| `Lao85` | ✅ | ✅ | ✅ |\n",
+    "| `TensionSpline` | ✅ | ❌ | ❌ |\n",
+    "| `GeneralPprimeFFprime` | ✅ | ❌ | ❌ |\n",
+    "\n",
+    "In this notebook, we will make use of the `ConstrainPaxisIp` (and `ConstrainBetapIp`) profiles (see [Jeon (2015)](https://link.springer.com/article/10.3938/jkps.67.843)). Others will be utilised in later notebooks. If there is a profile parameterisation you require that does not exist, please do create an issue.\n",
     "\n",
     "Both `ConstrainPaxisIp` and `ConstrainBetapIp` are parameterised as follows:\n",
     " $$J_{p}(\\psi, R, Z) = \\lambda\\big[ \\beta_{0} \\frac{R}{R_{0}} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} + (1-\\beta_{0}) \\frac{R_0}{R} \\left( 1-\\tilde{\\psi}^{\\alpha_m} \\right)^{\\alpha_n} \\big] \\quad (R,Z) \\in \\Omega_p,$$\n",

examples/example02 - static_forward_solve_MASTU.ipynb

@@ -599,6 +599,117 @@
     "print(f\"Original beta_0 = {profiles_topeol.Beta0} vs. Fitted from Lao85 = {beta_0}.\")"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Supposing an explicit parameterisation of the profiles is unknown and you have access to data describing the profiles. The `GeneralPprimeFFprime` profiles class enables you to use the raw data. \n",
+    "\n",
+    "As always the plasma current density is given by:\n",
+    "\n",
+    " $$J_{p}(\\psi, R, Z) = \\lambda\\big[ \\frac{R}{R_{0}} p'(\\tilde{\\psi}) + \\frac{R_0}{R} \\frac{1}{\\mu_0} F F'(\\tilde{\\psi}) \\big] \\quad (R,Z) \\in \\Omega_p, $$\n",
+    "\n",
+    "where the pressure, $p(\\tilde{\\psi})$, and toroidal magnetic field, $F(\\tilde{\\psi})$, profiles are now given by data arrays. \n",
+    "\n",
+    "By passing these data arrays either directly as `pprime_data` or `p_data` (and also `ffprime_data` or `f_data`), FreeGSNKE will interpolate the data and use them directly in the plasma current density above. \n",
+    "\n",
+    "The required parameters are:\n",
+    "- `Ip` (total plasma current).\n",
+    "- `fvac` ($rB_{tor}$,  vacuum toroidal field strength).\n",
+    "- `psi_n` ($\\tilde{\\psi}$, normalised poloidal flux at which data arrays corresponding, should be increasing).\n",
+    "- `pprime_data` or `p_data` (array of $p'(\\tilde{\\psi})$ or $p(\\tilde{\\psi})$ values at the normalised poloidal flux values in `psi_n`).\n",
+    "- `ffprime_data` or `f_data` (array of $FF'(\\tilde{\\psi})$ or $F(\\tilde{\\psi})$ values at the normalised poloidal flux values in `psi_n`).\n",
+    "- `Ip_logic` (if False, `Ip` is not used, if True, `Ip` is used to normalise $J_p$ and find $\\lambda$)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from freegsnke.jtor_update import GeneralPprimeFFprime\n",
+    "\n",
+    "psi_n = np.linspace(0,1,11)\n",
+    "\n",
+    "pprime_data = np.array([208061.72441505, 204110.95531309, 194359.61041113, 179811.09692964,\n",
+    "       161053.35395914, 138602.79928382, 113000.64896365,  84893.04354249,\n",
+    "        55164.04720078,  25298.53358472,      0.        ])\n",
+    "\n",
+    "ffprime_data = np.array([0.58961295, 0.57841712, 0.5507834 , 0.50955529, 0.45639891,\n",
+    "       0.3927777 , 0.32022539, 0.24057302, 0.1563259 , 0.07169191,\n",
+    "       0.        ])\n",
+    "\n",
+    "p_data = np.array([126175.75656716, 105510.77120406,  85544.30596718,  66798.45790761,\n",
+    "        49722.32952766,  34710.90959473,  22106.93373   ,  12194.53794856,\n",
+    "         5182.93833357,   1168.13190975,      0.        ])\n",
+    "\n",
+    "f_data = np.array([0.98240673, 0.92086948, 0.85722645, 0.79283784, 0.72925299,\n",
+    "       0.66837871, 0.61261316, 0.56490234, 0.52855955, 0.50657733,\n",
+    "       0.5       ])\n",
+    "\n",
+    "\n",
+    "# the following three setups are equivalent\n",
+    "profiles_general = GeneralPprimeFFprime(\n",
+    "    eq=eq,\n",
+    "    Ip=6e5,\n",
+    "    fvac=0.5,\n",
+    "    psi_n=psi_n,\n",
+    "    pprime_data=pprime_data,\n",
+    "    ffprime_data=ffprime_data,\n",
+    "    p_data=p_data,\n",
+    "    f_data=f_data,\n",
+    ")\n",
+    "\n",
+    "# profiles_general = GeneralPprimeFFprime(\n",
+    "#     eq=eq,\n",
+    "#     Ip=6e5,\n",
+    "#     fvac=0.5,\n",
+    "#     psi_n=psi_n,\n",
+    "#     pprime_data=pprime_data,\n",
+    "#     ffprime_data=ffprime_data,\n",
+    "#     p_data=None,\n",
+    "#     f_data=None,\n",
+    "# )\n",
+    "\n",
+    "# profiles_general = GeneralPprimeFFprime(\n",
+    "#     eq=eq,\n",
+    "#     Ip=6e5,\n",
+    "#     fvac=0.5,\n",
+    "#     psi_n=psi_n,\n",
+    "#     pprime_data=None,\n",
+    "#     ffprime_data=None,\n",
+    "#     p_data=p_data,\n",
+    "#     f_data=f_data,\n",
+    "# )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# instatiate new equilibrium object\n",
+    "eq_general = deepcopy(eq)\n",
+    "\n",
+    "# call solver with new profile object\n",
+    "GSStaticSolver.solve(eq=eq_general, \n",
+    "                     profiles=profiles_general, \n",
+    "                     constrain=None, \n",
+    "                     target_relative_tolerance=1e-9)\n",
+    "\n",
+    "\n",
+    "# plot the resulting equilbria \n",
+    "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n",
+    "ax1.grid(True, which='both')\n",
+    "eq_general.plot(axis=ax1, show=False)\n",
+    "eq_general.tokamak.plot(axis=ax1, show=False)\n",
+    "ax1.set_xlim(0.1, 2.15)\n",
+    "ax1.set_ylim(-2.25, 2.25)\n",
+    "plt.tight_layout()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

examples/example03 - extracting_equilibrium_quantites.ipynb

@@ -906,7 +906,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot the input p' and FF' profiles\n",
+    "# plot the p' and FF' profiles\n",
     "\n",
     "psi_n = eq.psiN_1D(N=65)\n",
     "\n",
@@ -929,15 +929,21 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot q profile\n",
+    "# plot the p and F profiles\n",
     "\n",
-    "psi_n = np.linspace(0.01,0.99,65)  # values of q at 0 and 1 can be problematic\n",
+    "psi_n = eq.psiN_1D(N=65)\n",
     "\n",
-    "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n",
+    "fig1, (ax1, ax2) = plt.subplots(1, 2, figsize=(15,6), dpi=80)\n",
     "ax1.grid(zorder=0, alpha=0.75)\n",
-    "ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax1.plot(psi_n, profiles.pressure(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
     "ax1.set_xlabel(r'$\\hat{\\psi}$')\n",
-    "ax1.set_ylabel(r\"$q(\\hat{\\psi})$\")\n"
+    "ax1.set_ylabel(r\"$p(\\hat{\\psi})$\")\n",
+    "ax1.ticklabel_format(axis='y', scilimits=(0,0))\n",
+    "\n",
+    "ax2.grid(zorder=0, alpha=0.75)\n",
+    "ax2.plot(psi_n, profiles.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax2.set_xlabel(r'$\\hat{\\psi}$')\n",
+    "ax2.set_ylabel(r\"$F(\\hat{\\psi})$\")\n"
    ]
   },
   {
@@ -946,14 +952,15 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# plot fpol\n",
+    "# plot q profile\n",
+    "\n",
+    "psi_n = np.linspace(0.01,0.99,65)  # values of q at 0 and 1 can be problematic\n",
     "\n",
     "fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n",
     "ax1.grid(zorder=0, alpha=0.75)\n",
-    "ax1.plot(psi_n, eq.fpol(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax1.plot(psi_n, eq.q(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
     "ax1.set_xlabel(r'$\\hat{\\psi}$')\n",
-    "ax1.set_ylabel(r\"$fpol(\\hat{\\psi})$\")\n",
-    "\n"
+    "ax1.set_ylabel(r\"$q(\\hat{\\psi})$\")\n"
    ]
   },
   {