Merge pull request #30 from FusionComputingLab/flux_averaging_example

Added example of flux averaging a 2D field

Author: nicamo

.github/workflows/notebooks.yml

@@ -8,6 +8,9 @@ on:
 
 jobs:
   notebooks:
+    strategy:
+      fail-fast: false
+    timeout-minutes: 60
     if: ${{ github.event.label.name == 'ready-for-final-tests' }}
     runs-on: ubuntu-latest
     steps:

docs/conf.py

@@ -130,4 +130,4 @@
     "dollarmath",
 ]
 
-nb_execution_timeout = 600
+nb_execution_timeout = 1200

examples/example03 - extracting_equilibrium_quantites.ipynb

@@ -69,6 +69,7 @@
     "\n",
     "for key in current_values.keys():\n",
     "    eq.tokamak[key].current = current_values[key]\n",
+    "eq.tokamak[\"P6\"].current += 100\n",
     "\n",
     "# carry out forward solve\n",
     "GSStaticSolver.solve(eq=eq, \n",
@@ -920,7 +921,30 @@
     "ax2.grid(zorder=0, alpha=0.75)\n",
     "ax2.plot(psi_n, profiles.ffprime(psi_n), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
     "ax2.set_xlabel(r'$\\hat{\\psi}$')\n",
-    "ax2.set_ylabel(r\"$FF'(\\hat{\\psi})$\")\n"
+    "ax2.set_ylabel(r\"$FF'(\\hat{\\psi})$\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot the p and F profiles\n",
+    "\n",
+    "psi_n = eq.psiN_1D(N=65)\n",
+    "\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, 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\"$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"
    ]
   },
   {
@@ -996,20 +1020,6 @@
     "plt.tight_layout()"
    ]
   },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "# # plot 1D_jtor\n",
-    "# fig1, ax1 = plt.subplots(1, 1, figsize=(6,6), dpi=80)\n",
-    "# ax1.grid(zorder=0, alpha=0.75)\n",
-    "# ax1.plot(eq.psiN_1D(N=101), eq.jtor_1D(N=101), color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
-    "# ax1.set_xlabel(r'$\\hat{\\psi}$')\n",
-    "# ax1.set_ylabel(r\"$J_{tor}(\\hat{\\psi})$\")"
-   ]
-  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -1054,6 +1064,96 @@
     "\n",
     "plt.tight_layout()"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Flux averaged quantities\n",
+    "\n",
+    "The `equilibrium` class provides a method to calculate the \"flux averaged\" value of a user-defined 2D scalar field $f(R,Z)$, using [line integrals](https://tutorial.math.lamar.edu/classes/calciii/LineIntegralsPtI.aspx), on a given (normalised) flux surface of $\\psi_n$ (within the last closed flux surface). The flux average $\\langle f \\rangle$ is given by\n",
+    "\n",
+    "$$\n",
+    "\\langle f \\rangle (\\psi_n) = \\frac{ \\int_{C(\\psi_n)} \\frac{f(R,Z)}{B_{\\text{pol}}(R,Z)}\\, ds}{ \\int_{C(\\psi_n)} \\frac{1}{B_{\\text{pol}}(R,Z)} \\, ds },\n",
+    "$$\n",
+    "\n",
+    "where:\n",
+    "- $f(R,Z)$ = 2D scalar field function (e.g. the plasma current density function $J_p(R,Z)$).\n",
+    "- $\\psi_n$ = value of normalised flux at which to evaluate line integrals.\n",
+    "- $C(\\psi_n)$ = curve of $(R, Z)$ points satisfying $\\psi_n = \\text{const}$ (i.e. a flux contour).\n",
+    "- $B_{\\text{pol}}(R,Z)$ = 2D scalar poloidal magnetic field function.\n",
+    "- $ds$ = the arc length element (where $ds = \\sqrt{(dR/dl)^2 + (dZ/dl)^2} dl$ and $l \\in [0,L]$ is a parameterised length going from the beginning to the end of the contour).\n",
+    "\n",
+    "The definition of the flux average was taken from [Song et al. (2024)](https://www.mdpi.com/2571-6182/7/4/45)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Here we will calculate the flux averaged values of the plasma current density by first defining a function that returns the current density at arbitrary $(R,Z)$ locations. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "from scipy.interpolate import RectBivariateSpline\n",
+    "\n",
+    "def f(R,Z):\n",
+    "    jtor = RectBivariateSpline(eq.R_1D, eq.Z_1D, eq._profiles.jtor)\n",
+    "    return jtor(R, Z, grid=False)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Next we can call the method in the equilibrium object and plot the results. Given the notation above, we note that $\\psi_n$ and $\\hat{\\psi}$ are equivalent. "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# call the method\n",
+    "flux_averaged_jtor, psi_n = eq.flux_averaged_function(\n",
+    "    f=f,\n",
+    "    psi_n=np.linspace(0.0,1.0,101)\n",
+    "    )"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot\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, flux_averaged_jtor, color='k', linewidth=1, marker='x', markersize=2, zorder=10)\n",
+    "ax1.set_xlabel(r'$\\hat{\\psi}$')\n",
+    "ax1.set_ylabel(r\"$\\langle J_p \\rangle (\\hat{\\psi})$\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# # for example you could flux average other quantities of interest\n",
+    "\n",
+    "# # 1/R\n",
+    "# def f(R,Z):\n",
+    "#     g = RectBivariateSpline(eq.R_1D, eq.Z_1D, 1/eq.R)\n",
+    "#     return g(R, Z, grid=False)"
+   ]
   }
  ],
  "metadata": {

requirements-freegs4e.txt

@@ -1 +1 @@
-freegs4e~=0.10
+freegs4e~=0.11