From 2a579c308b318a7c5a5372209dbe1e005eb86458 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 26 Mar 2026 15:54:55 +0000 Subject: [PATCH 01/11] new example using the Anamak tokamak from FGE paper --- examples/example07a - Anamak.ipynb | 786 +++++++++++++++++++++++++++++ 1 file changed, 786 insertions(+) create mode 100644 examples/example07a - Anamak.ipynb diff --git a/examples/example07a - Anamak.ipynb b/examples/example07a - Anamak.ipynb new file mode 100644 index 00000000..028e1535 --- /dev/null +++ b/examples/example07a - Anamak.ipynb @@ -0,0 +1,786 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example: \"Anamak\" tokamak\n", + "---\n", + "\n", + "Here, we generate a few different equilibria using the \"Anamak\" example tokamak from the [FGE code paper](https://arxiv.org/abs/2512.06847).\n", + "\n", + "The equilbirium\\profile parameters are **completely made up** based on some equilbiria in the above paper - please experiment on your own and change them to more realistic values as you please!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create the machine object\n", + "\n", + "The Anamak machine is a useful toy tokamak because the geometry can be varied in a simple and controlled way. In particular, the number and placement of active coils, passive structures, and the vessel can all be adjusted using a small set of geometric parameters.\n", + "\n", + "Following the parameterisation used in the reference above, the $N_{\\mathrm{coils}}$ **active coils** are placed on an **ellipse** centred at $(R_0, 0)$ in the poloidal $(R, Z)$ plane. Their coordinates are defined by\n", + "\n", + "$$\n", + "R_{\\mathrm{coil}}^{(i)} = R_0 + \\hat{R}_{\\mathrm{coil}} \\cos(\\theta_i),\n", + "$$\n", + "\n", + "$$\n", + "Z_{\\mathrm{coil}}^{(i)} = \\kappa \\, \\hat{R}_{\\mathrm{coil}} \\sin(\\theta_i),\n", + "$$\n", + "\n", + "where\n", + "- $R_0$ is the major radius at which the ellipse is centred.\n", + "- $\\hat{R}_{\\mathrm{coil}}$ is the minor radius of the coils. \n", + "- $\\bold{\\theta}$ is a vector (length $N_{\\mathrm{coils}}$) of poloidal angles $\\theta \\in [0, 2\\pi)$. \n", + "- $\\kappa$ is the \"elongation\" of the ellipse.\n", + "\n", + "Similarly, the $N_{\\mathrm{passives}}$ **passive structures** are placed on an ellipse at the same centre:\n", + "\n", + "$$\n", + "R_{\\mathrm{passives}}^{(i)} = R_0 + \\hat{R}_{\\mathrm{passive}} \\cos(\\theta_i),\n", + "$$\n", + "\n", + "$$\n", + "Z_{\\mathrm{passives}}^{(i)} = \\kappa \\, \\hat{R}_{\\mathrm{passive}} \\sin(\\theta_i),\n", + "$$\n", + "\n", + "but with a potentially different minor radius $ \\hat{R}_{\\mathrm{passive}}$. \n", + "\n", + "Finally the $N_{\\mathrm{limiter}}$ limiter/wall points are set up similarly:\n", + "\n", + "$$\n", + "R_{\\mathrm{limiter}}^{(i)} = R_0 + \\hat{R}_{\\mathrm{limiter}} \\cos(\\theta_i),\n", + "$$\n", + "\n", + "$$\n", + "Z_{\\mathrm{limiter}}^{(i)} = \\kappa \\, \\hat{R}_{\\mathrm{limiter}} \\sin(\\theta_i).\n", + "$$\n", + "\n", + "We should note that given we are only doing static solves below, the passive structures are not really used here." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In the next cell, we can choose the parameters for building the Anamak." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "# major radius of the machine\n", + "R0 = 1.0\n", + "\n", + "# \"elongation\" of the machine (same for actives, passives, and limiter here)\n", + "kappa = 1.0 \n", + "\n", + "# for each of the coils, passives, and limiter, define:\n", + "# -> N = the number of each you desire\n", + "# -> R = the minor radius of its ellipse\n", + "\n", + "# actives parameters\n", + "N_actives = 8\n", + "R_coil = 0.7\n", + "\n", + "# passives parameters\n", + "N_passives = 50\n", + "R_passive = 0.5\n", + "\n", + "# limiter parameters\n", + "N_limiter = 101\n", + "R_limiter = 0.45" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The following functions are used to build the Anamak." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# builds an ellipse of (R,Z) points\n", + "def build_ellipse(\n", + " R0, # major radius\n", + " N, # number of points in [0,2\\pi)\n", + " R_minor, # minor radius\n", + " kappa, # \"elongation\"\n", + " ):\n", + " \n", + " # poloidal locations\n", + " theta = np.arange(0, 2*np.pi, 2*np.pi/N)\n", + " \n", + " # locations\n", + " r_locs = R0 + R_minor*np.cos(theta)\n", + " z_locs = R_minor*kappa*np.sin(theta)\n", + "\n", + " return r_locs, z_locs\n", + "\n", + "# builds active coils using (R,Z) locations of filaments\n", + "def build_actives(\n", + " r_locs, \n", + " z_locs, \n", + " ):\n", + "\n", + " # populate active coils dictionary with single filament coils (point-sources)\n", + " active_coils_data = {}\n", + " for i, (r, z) in enumerate(zip(r_locs, z_locs)):\n", + " active_coils_data[i] = {\n", + " \"R\": [r],\n", + " \"Z\": [z],\n", + " \"dR\": 0.05, # size only used for plotting in the simulations below\n", + " \"dZ\": 0.05, # size only used for plotting in the simulations below\n", + " \"resistivity\": 1.55e-8, # not used below\n", + " \"polarity\": 1, # default\n", + " \"multiplier\": 1, # single turn coils\n", + " }\n", + " \n", + " return active_coils_data\n", + "\n", + "# builds passive structures using (R,Z) locations of filaments\n", + "def build_passives(\n", + " r_locs, \n", + " z_locs, \n", + " ):\n", + "\n", + " # populate passive coils list with single filament passives (point-sources)\n", + " passive_coils_data = []\n", + " for r, z in zip(r_locs, z_locs):\n", + " passive_coils_data.append({\n", + " \"R\": r,\n", + " \"Z\": z,\n", + " \"dR\": 0.02, # size only used for plotting in the simulations below\n", + " \"dZ\": 0.02, # size only used for plotting in the simulations below\n", + " \"resistivity\": 5.5e-7 # resistivity of material (steel)\n", + " })\n", + " \n", + " return passive_coils_data\n", + "\n", + "# builds limiter geometry using (R,Z) locations\n", + "def build_limiter(\n", + " r_locs, \n", + " z_locs, \n", + " ):\n", + "\n", + " # populate limiter list\n", + " limiter_data = []\n", + " for r, z in zip(r_locs, z_locs):\n", + " limiter_data.append({\"R\": r, \"Z\": z})\n", + " \n", + " return limiter_data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build the active coils data\n", + "r_active_locs, z_active_locs = build_ellipse(\n", + " R0=R0, \n", + " N=N_actives, \n", + " R_minor=1.4*R_passive,\n", + " kappa=kappa,\n", + ")\n", + "\n", + "active_coils_data = build_actives(\n", + " r_locs=r_active_locs, \n", + " z_locs=z_active_locs\n", + " )\n", + "\n", + "# build the passive structures data\n", + "r_passive_locs, z_passive_locs = build_ellipse(\n", + " R0=R0, \n", + " N=N_passives, \n", + " R_minor=R_passive,\n", + " kappa=kappa,\n", + ")\n", + "\n", + "passive_coils_data = build_passives(\n", + " r_locs=r_passive_locs, \n", + " z_locs=z_passive_locs\n", + " )\n", + "\n", + "# build the limiter data\n", + "r_limiter_locs, z_limiter_locs = build_ellipse(\n", + " R0=R0, \n", + " N=N_limiter, \n", + " R_minor=R_limiter,\n", + " kappa=kappa,\n", + ")\n", + "\n", + "limiter_data = build_limiter(\n", + " r_locs=r_limiter_locs, \n", + " z_locs=z_limiter_locs\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build the Anamak\n", + "from freegsnke import build_machine\n", + "tokamak = build_machine.tokamak(\n", + " active_coils_data=active_coils_data,\n", + " passive_coils_data=passive_coils_data,\n", + " limiter_data=limiter_data,\n", + " wall_data=limiter_data,\n", + ")\n", + "\n", + "# plot the machine\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "tokamak.plot(axis=ax1,show=False) \n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) \n", + "ax1.set_xlabel(\"Major radius [m]\")\n", + "ax1.set_ylabel(\"Height [m]\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate an equilibrium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import equilibrium_update\n", + "\n", + "eq = equilibrium_update.Equilibrium(\n", + " tokamak=tokamak, # provide tokamak object\n", + " Rmin=0.1, Rmax=1.9, # radial range\n", + " Zmin=-0.8, Zmax=0.8, # vertical range\n", + " nx=65, # 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", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate a profile object\n", + "\n", + "We use the Lao profile object here, providing constant pressure and toroidal field profiles for simplicity. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke.jtor_update import Lao85\n", + "profiles = Lao85(\n", + " eq=eq, # eq object\n", + " Ip=2e5, # plasma current\n", + " fvac=0.5*R0, # toroidal field * major radius\n", + " alpha=[58213.6], # linear profile\n", + " beta=[0.582136] # linear profile\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the static nonlinear solver" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import GSstaticsolver\n", + "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Inverse solver constraints\n", + "\n", + "In this first example, we'll set up some simple constraints to generate an up/down symmetric double null plasma (with relatively high triangularity). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# null points constraints\n", + "Rx = 0.8 # X-point radius\n", + "Zx = 0.35 # X-point height\n", + "null_points = [[Rx, Rx], [Zx, -Zx]]\n", + "\n", + "# isoflux constraints (including null points here)\n", + "Rout = 1.22 # outboard midplane radius\n", + "isoflux_set = np.array([\n", + " [\n", + " [Rx, Rx, Rout], \n", + " [Zx, -Zx, 0.0]\n", + " ]]\n", + " )\n", + " \n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(\n", + " null_points=null_points,\n", + " isoflux_set=isoflux_set,\n", + " weight_isoflux=1.0,\n", + " weight_nulls=0.8, # slightly lower weight on null points to get correct shape\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The inverse solve\n", + "\n", + "Solve using the inverse solver with forced up/down symmetry. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-5,\n", + " target_relative_psit_update=1e-3,\n", + " verbose=True,\n", + " l2_reg=np.array([1e-12]*N_actives),\n", + " force_up_down_symmetric=True,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# refine with a forward solve (now the currents are known)\n", + "GSStaticSolver.solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=None, \n", + " target_relative_tolerance=1e-9,\n", + " verbose=False, # print output\n", + " )" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Using only a few constraints, we can see that the inverse solver has found the desired plasma. It is almost perfectly up/down symmetric and has high triangularity. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Triangularity = {eq.triangularity()}\")\n", + "print(f\"Avg. Z (Jtor) position = {eq.Zcurrent()} [m]\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the equilibrium\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next, we'll use the same machine description to make a negtive triangularity plasma (again up/down symmetric). " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# null points constraints\n", + "Rx = 1.1 # X-point radius\n", + "Zx = 0.35 # X-point height\n", + "null_points = [[Rx, Rx], [Zx, -Zx]]\n", + "\n", + "# isoflux constraints (including null points here)\n", + "Rin = 0.8 # inboard midplane radius\n", + "isoflux_set = np.array([\n", + " [\n", + " [Rx, Rx, Rin], \n", + " [Zx, -Zx, 0.0]\n", + " ]]\n", + " )\n", + " \n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(\n", + " null_points=null_points,\n", + " isoflux_set=isoflux_set,\n", + " weight_isoflux=1.0,\n", + " weight_nulls=0.8, # slightly lower weight on null points to get correct shape\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# solve\n", + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-5,\n", + " target_relative_psit_update=1e-3,\n", + " verbose=True,\n", + " l2_reg=np.array([1e-14]*N_actives),\n", + " force_up_down_symmetric=True,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# refine with a forward solve (now the currents are known)\n", + "GSStaticSolver.solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=None, \n", + " target_relative_tolerance=1e-9,\n", + " verbose=False, # print output\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the equilibrium\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now lets look at changing the shape of the Anamak.\n", + "\n", + "We'll now increase the number of active coils and increase the \"elongation\" of the device. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# \"elongation\" of the machine (same for actives, passives, and limiter here)\n", + "kappa = 1.8\n", + "\n", + "# for each of the coils, passives, and limiter, define:\n", + "# -> N = the number of each you desire\n", + "# -> R = the minor radius of its ellipse\n", + "\n", + "# actives parameters\n", + "N_actives = 24\n", + "R_coil = 0.7" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build the active coils data\n", + "r_active_locs, z_active_locs = build_ellipse(\n", + " R0=R0, \n", + " N=N_actives, \n", + " R_minor=1.4*R_passive,\n", + " kappa=kappa,\n", + ")\n", + "\n", + "active_coils_data = build_actives(\n", + " r_locs=r_active_locs, \n", + " z_locs=z_active_locs\n", + " )\n", + "\n", + "# build the passive structures data\n", + "r_passive_locs, z_passive_locs = build_ellipse(\n", + " R0=R0, \n", + " N=N_passives, \n", + " R_minor=R_passive,\n", + " kappa=kappa,\n", + ")\n", + "\n", + "passive_coils_data = build_passives(\n", + " r_locs=r_passive_locs, \n", + " z_locs=z_passive_locs\n", + " )\n", + "\n", + "# build the limiter data\n", + "r_limiter_locs, z_limiter_locs = build_ellipse(\n", + " R0=R0, \n", + " N=N_limiter, \n", + " R_minor=R_limiter,\n", + " kappa=kappa,\n", + ")\n", + "\n", + "limiter_data = build_limiter(\n", + " r_locs=r_limiter_locs, \n", + " z_locs=z_limiter_locs\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build the Anamak\n", + "from freegsnke import build_machine\n", + "tokamak = build_machine.tokamak(\n", + " active_coils_data=active_coils_data,\n", + " passive_coils_data=passive_coils_data,\n", + " limiter_data=limiter_data,\n", + " wall_data=limiter_data,\n", + ")\n", + "\n", + "# plot the machine\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "tokamak.plot(axis=ax1,show=False) \n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) \n", + "ax1.set_xlabel(\"Major radius [m]\")\n", + "ax1.set_ylabel(\"Height [m]\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Re-initialise the equilibrium and profiles objects: here we've increased the vertical size of the computational grid and increased the plasma current. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import equilibrium_update\n", + "\n", + "eq = equilibrium_update.Equilibrium(\n", + " tokamak=tokamak, # provide tokamak object\n", + " Rmin=0.1, Rmax=1.9, # radial range\n", + " Zmin=-1.5, Zmax=1.5, # vertical range\n", + " nx=65, # 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", + "from freegsnke.jtor_update import Lao85\n", + "profiles = Lao85(\n", + " eq=eq, # eq object\n", + " Ip=5e5, # plasma current\n", + " fvac=0.5*R0, # toroidal field * major radius\n", + " alpha=[58213.6], # linear profile\n", + " beta=[0.582136] # linear profile\n", + ")\n", + "\n", + "from freegsnke import GSstaticsolver\n", + "GSStaticSolver = GSstaticsolver.NKGSsolver(eq)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's try to create a simple up/down symmetric double null plasma in this new machine. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# null points constraints\n", + "Rx = 0.8 # X-point radius\n", + "Zx = 0.65 # X-point height\n", + "null_points = [[Rx, Rx], [Zx, -Zx]]\n", + "\n", + "# isoflux constraints (including null points here)\n", + "Rout = 1.3 # outboard midplane radius\n", + "isoflux_set = np.array([\n", + " [\n", + " [Rx, Rx, Rout], \n", + " [Zx, -Zx, 0.0]\n", + " ]]\n", + " )\n", + " \n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(\n", + " null_points=null_points,\n", + " isoflux_set=isoflux_set,\n", + " weight_isoflux=1.0,\n", + " weight_nulls=0.7, # slightly lower weight on null points to get correct shape\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# solve\n", + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-5,\n", + " target_relative_psit_update=1e-3,\n", + " verbose=True,\n", + " l2_reg=np.array([1e-12]*N_actives),\n", + " force_up_down_symmetric=True,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# refine with a forward solve (now the currents are known)\n", + "GSStaticSolver.solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=None, \n", + " target_relative_tolerance=1e-9,\n", + " verbose=False, # print output\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the equilibrium\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints" + ] + }, + { + "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": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "freegsnke4e", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.16" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 312bf0064131df3e933d719dc6ee058008c7f264 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 26 Mar 2026 15:55:18 +0000 Subject: [PATCH 02/11] re-naming SPARC example --- examples/example07b - SPARC.ipynb | 261 ++++++++++++++++++++++++++++++ 1 file changed, 261 insertions(+) create mode 100644 examples/example07b - SPARC.ipynb diff --git a/examples/example07b - SPARC.ipynb b/examples/example07b - SPARC.ipynb new file mode 100644 index 00000000..ba8360e5 --- /dev/null +++ b/examples/example07b - SPARC.ipynb @@ -0,0 +1,261 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example: Static inverse free-boundary equilibrium calculations (in SPARC)\n", + "\n", + "---\n", + "\n", + "Here we will generate an equilibrium (find coil currents with the inverse solver) in a SPARC-like tokamak. \n", + "\n", + "The machine description comes from files located [here](https://github.com/cfs-energy/SPARCPublic).\n", + "\n", + "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create the machine object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build machine\n", + "from freegsnke import build_machine\n", + "tokamak = build_machine.tokamak(\n", + " active_coils_path=f\"../machine_configs/SPARC/SPARC_active_coils.pickle\",\n", + " passive_coils_path=f\"../machine_configs/SPARC/SPARC_passive_coils.pickle\",\n", + " limiter_path=f\"../machine_configs/SPARC/SPARC_limiter.pickle\",\n", + " wall_path=f\"../machine_configs/SPARC/SPARC_wall.pickle\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate an equilibrium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import equilibrium_update\n", + "\n", + "eq = equilibrium_update.Equilibrium(\n", + " tokamak=tokamak, # provide tokamak object\n", + " Rmin=1.1, Rmax=2.7, # radial range\n", + " Zmin=-1.8, Zmax=1.8, # 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", + " # psi=plasma_psi\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate a profile object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialise the profiles\n", + "from freegsnke.jtor_update import ConstrainPaxisIp\n", + "profiles = ConstrainPaxisIp(\n", + " eq=eq, # equilibrium object\n", + " paxis=5e4, # pressure on axis\n", + " Ip=8.7e6, # 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", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the static nonlinear solver" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import GSstaticsolver\n", + "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Constraints" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Rx = 1.55 # X-point radius\n", + "Zx = 1.15 # X-point height\n", + "\n", + "# set desired null_points locations\n", + "# this can include X-point and O-point locations\n", + "null_points = [[Rx, Rx], [Zx, -Zx]]\n", + "\n", + "Rout = 2.4 # outboard midplane radius\n", + "Rin = 1.3 # inboard midplane radius\n", + "\n", + "# set desired isoflux constraints with format \n", + "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", + "# with each isoflux_i = [R_coords, Z_coords]\n", + "isoflux_set = np.array([[[Rx, Rx, Rout, Rin, 1.7, 1.7], [Zx, -Zx, 0.0, 0.0, 1.5, -1.5]]])\n", + " \n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(null_points=null_points,\n", + " isoflux_set=isoflux_set)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The inverse solve" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-6,\n", + " target_relative_psit_update=1e-3,\n", + " verbose=False, # print output\n", + " l2_reg=np.array([1e-16]*10 + [1e-5]), \n", + " )\n", + " \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# refine with a forward solve (now the currents are known)\n", + "GSStaticSolver.solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=None, \n", + " target_relative_tolerance=1e-9,\n", + " verbose=False, # print output\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints\n", + "ax1.set_xlim(1.0, 3.0)\n", + "ax1.set_ylim(-2.0, 2.0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eq.tokamak.getCurrents()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "freegsnke4e", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.16" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 7ee2c9e95888773b1852bcd1f07fac6645169bed Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 26 Mar 2026 15:56:44 +0000 Subject: [PATCH 03/11] removing old SPARC example --- ...ample07 - static_inverse_solve_SPARC.ipynb | 261 ------------------ 1 file changed, 261 deletions(-) delete mode 100644 examples/example07 - static_inverse_solve_SPARC.ipynb diff --git a/examples/example07 - static_inverse_solve_SPARC.ipynb b/examples/example07 - static_inverse_solve_SPARC.ipynb deleted file mode 100644 index ba8360e5..00000000 --- a/examples/example07 - static_inverse_solve_SPARC.ipynb +++ /dev/null @@ -1,261 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Example: Static inverse free-boundary equilibrium calculations (in SPARC)\n", - "\n", - "---\n", - "\n", - "Here we will generate an equilibrium (find coil currents with the inverse solver) in a SPARC-like tokamak. \n", - "\n", - "The machine description comes from files located [here](https://github.com/cfs-energy/SPARCPublic).\n", - "\n", - "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Import packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numpy as np" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Create the machine object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# build machine\n", - "from freegsnke import build_machine\n", - "tokamak = build_machine.tokamak(\n", - " active_coils_path=f\"../machine_configs/SPARC/SPARC_active_coils.pickle\",\n", - " passive_coils_path=f\"../machine_configs/SPARC/SPARC_passive_coils.pickle\",\n", - " limiter_path=f\"../machine_configs/SPARC/SPARC_limiter.pickle\",\n", - " wall_path=f\"../machine_configs/SPARC/SPARC_wall.pickle\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", - "\n", - "ax1.grid(zorder=0, alpha=0.75)\n", - "ax1.set_aspect('equal')\n", - "tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", - "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Instantiate an equilibrium" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke import equilibrium_update\n", - "\n", - "eq = equilibrium_update.Equilibrium(\n", - " tokamak=tokamak, # provide tokamak object\n", - " Rmin=1.1, Rmax=2.7, # radial range\n", - " Zmin=-1.8, Zmax=1.8, # 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", - " # psi=plasma_psi\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Instantiate a profile object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# initialise the profiles\n", - "from freegsnke.jtor_update import ConstrainPaxisIp\n", - "profiles = ConstrainPaxisIp(\n", - " eq=eq, # equilibrium object\n", - " paxis=5e4, # pressure on axis\n", - " Ip=8.7e6, # 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", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Load the static nonlinear solver" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke import GSstaticsolver\n", - "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Constraints" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Rx = 1.55 # X-point radius\n", - "Zx = 1.15 # X-point height\n", - "\n", - "# set desired null_points locations\n", - "# this can include X-point and O-point locations\n", - "null_points = [[Rx, Rx], [Zx, -Zx]]\n", - "\n", - "Rout = 2.4 # outboard midplane radius\n", - "Rin = 1.3 # inboard midplane radius\n", - "\n", - "# set desired isoflux constraints with format \n", - "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", - "# with each isoflux_i = [R_coords, Z_coords]\n", - "isoflux_set = np.array([[[Rx, Rx, Rout, Rin, 1.7, 1.7], [Zx, -Zx, 0.0, 0.0, 1.5, -1.5]]])\n", - " \n", - "# instantiate the freegsnke constrain object\n", - "from freegsnke.inverse import Inverse_optimizer\n", - "constrain = Inverse_optimizer(null_points=null_points,\n", - " isoflux_set=isoflux_set)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### The inverse solve" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "GSStaticSolver.inverse_solve(eq=eq, \n", - " profiles=profiles, \n", - " constrain=constrain, \n", - " target_relative_tolerance=1e-6,\n", - " target_relative_psit_update=1e-3,\n", - " verbose=False, # print output\n", - " l2_reg=np.array([1e-16]*10 + [1e-5]), \n", - " )\n", - " \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# refine with a forward solve (now the currents are known)\n", - "GSStaticSolver.solve(eq=eq, \n", - " profiles=profiles, \n", - " constrain=None, \n", - " target_relative_tolerance=1e-9,\n", - " verbose=False, # print output\n", - " )" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", - "\n", - "ax1.grid(zorder=0, alpha=0.75)\n", - "ax1.set_aspect('equal')\n", - "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", - "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", - "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", - "constrain.plot(axis=ax1, show=False) # plots the contraints\n", - "ax1.set_xlim(1.0, 3.0)\n", - "ax1.set_ylim(-2.0, 2.0)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "eq.tokamak.getCurrents()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "freegsnke4e", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.16" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From ab14722ba616da753ad3938fb60418d5cc70202c Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 26 Mar 2026 15:57:24 +0000 Subject: [PATCH 04/11] removing ITER example while we fix --- ...xample08 - static_inverse_solve_ITER.ipynb | 331 ------------------ 1 file changed, 331 deletions(-) delete mode 100644 examples/example08 - static_inverse_solve_ITER.ipynb diff --git a/examples/example08 - static_inverse_solve_ITER.ipynb b/examples/example08 - static_inverse_solve_ITER.ipynb deleted file mode 100644 index 1249af08..00000000 --- a/examples/example08 - static_inverse_solve_ITER.ipynb +++ /dev/null @@ -1,331 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Example: Static inverse free-boundary equilibrium calculations (in ITER)\n", - "\n", - "---\n", - "\n", - "Here we will generate an equilibrium (find coil currents with the inverse solver) in an ITER-like tokamak. \n", - "\n", - "The machine description comes from files located [here](https://github.com/ProjectTorreyPines/FUSE.jl).\n", - "\n", - "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Import packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numpy as np" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Create the machine object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# build machine\n", - "from freegsnke import build_machine\n", - "tokamak = build_machine.tokamak(\n", - " active_coils_path=f\"../machine_configs/ITER/ITER_active_coils.pickle\",\n", - " passive_coils_path=f\"../machine_configs/ITER/ITER_passive_coils.pickle\",\n", - " limiter_path=f\"../machine_configs/ITER/ITER_limiter.pickle\",\n", - " wall_path=f\"../machine_configs/ITER/ITER_wall.pickle\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot the machine\n", - "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", - "plt.tight_layout()\n", - "\n", - "tokamak.plot(axis=ax1, show=False)\n", - "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", - "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", - "\n", - "ax1.grid(alpha=0.5)\n", - "ax1.set_aspect('equal')\n", - "ax1.set_xlim(1.1, 12.4)\n", - "ax1.set_ylim(-8.3, 8.3)\n", - "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", - "ax1.set_ylabel(r'Height, $Z$ [m]')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Instantiate an equilibrium" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke import equilibrium_update\n", - "\n", - "eq = equilibrium_update.Equilibrium(\n", - " tokamak=tokamak, # provide tokamak object\n", - " Rmin=3.2, Rmax=8.8, # radial range\n", - " Zmin=-5, Zmax=5, # 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", - " # psi=plasma_psi\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Instantiate a profile object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# initialise the profiles\n", - "from freegsnke.jtor_update import ConstrainBetapIp\n", - "profiles = ConstrainBetapIp(\n", - " eq=eq, # equilibrium object\n", - " betap=0.15, # poloidal beta\n", - " Ip=11e6, # plasma current\n", - " fvac=0.5, # fvac = rB_{tor}\n", - " alpha_m=2.0, # profile function parameter\n", - " alpha_n=1.0 # profile function parameter\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Load the static nonlinear solver" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke import GSstaticsolver\n", - "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Constraints" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Rx = 5.02 # X-point radius\n", - "Zx = -3.23 # X-point height\n", - "\n", - "Ro = 6.34 # X-point radius\n", - "Zo = 0.66 # X-point height\n", - "\n", - "# set desired null_points locations\n", - "# this can include X-point and O-point locations\n", - "null_points = [[Rx, Ro], [Zx, Zo]]\n", - "\n", - "# set desired isoflux constraints with format \n", - "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", - "# with each isoflux_i = [R_coords, Z_coords]\n", - "isoflux_set = np.array([[\n", - " [4.25455147,\n", - " 4.1881875,\n", - " 4.2625,\n", - " 4.45683769,\n", - " 4.78746942,\n", - " 5.31835938,\n", - " 5.91875,\n", - " 6.44275166,\n", - " 6.92285156,\n", - " 7.35441013,\n", - " 7.73027344,\n", - " 8.03046875,\n", - " 8.19517497,\n", - " 8.05673414,\n", - " 7.75308013,\n", - " 7.37358451,\n", - " 6.94355469,\n", - " 6.47773438,\n", - " 5.99121094,\n", - " 5.49433594,\n", - " Rx,\n", - " 4.79042969,\n", - " 4.56269531,\n", - " 4.36038615], \n", - "[0.0,\n", - " 1.13554688,\n", - " 2.2565134,\n", - " 3.16757813,\n", - " 3.80507812,\n", - " 4.0499021,\n", - " 3.93089003,\n", - " 3.665625,\n", - " 3.31076604,\n", - " 2.86875,\n", - " 2.3199601,\n", - " 1.62832118,\n", - " 0.657421875,\n", - " -0.35859375,\n", - " -1.05585937,\n", - " -1.59375,\n", - " -2.03492727,\n", - " -2.41356403,\n", - " -2.75622016,\n", - " -3.08699278,\n", - " Zx,\n", - " -2.40548044,\n", - " -1.55982142,\n", - " -0.67734375]\n", - " ]])\n", - " \n", - "\n", - "# instantiate the freegsnke constrain object\n", - "from freegsnke.inverse import Inverse_optimizer\n", - "constrain = Inverse_optimizer(null_points=null_points,\n", - " isoflux_set=isoflux_set)\n", - "\n", - "\n", - "# remove from the coils available for control the radial field coil \n", - "# eq.tokamak['VS3'].control = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### The inverse solve" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "GSStaticSolver.inverse_solve(eq=eq, \n", - " profiles=profiles, \n", - " constrain=constrain, \n", - " target_relative_tolerance=1e-4,\n", - " target_relative_psit_update=1e-3,\n", - " # max_solving_iterations=150,\n", - " # max_iter_per_update=10,\n", - " # max_rel_update_size=0.075,\n", - " # damping_factor=.99,\n", - " # max_rel_psit=.05,\n", - " verbose=True, # print output\n", - " l2_reg=1e-14,\n", - " )\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", - "\n", - "ax1.grid(zorder=0, alpha=0.75)\n", - "ax1.set_aspect('equal')\n", - "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", - "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", - "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", - "constrain.plot(axis=ax1, show=False) # plots the contraints\n", - "ax1.set_xlim(1.1, 12.4)\n", - "ax1.set_ylim(-8.3, 8.3)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "eq.tokamak.getCurrents()\n", - "\n", - "# # save coil currents to file\n", - "# import pickle\n", - "# with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", - "# pickle.dump(obj=inverse_current_values, file=f)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] - } - ], - "metadata": { - "kernelspec": { - "display_name": "freegsnke4e", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.10.16" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From 2ad02b29ccc74d66b7500dbb6b9982986352ed05 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 26 Mar 2026 16:00:26 +0000 Subject: [PATCH 05/11] updating examples README --- examples/README.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/README.md b/examples/README.md index 670e3004..ae1c2437 100644 --- a/examples/README.md +++ b/examples/README.md @@ -12,8 +12,8 @@ These example Jupyter notebooks are intended to be the **first port of call for | Example 04 | Learn how to use the magnetic probes object. | Anyone | | Example 05 | Learn how to use the evolutive solver to simulate time-dependent equilibria. | Anyone | | Example 06a/b | Simulate (static) MAST-U equilibria over an entire shot using inputs from EFIT++ (requires internal UKAEA MAST-U database). | UKAEA employees + collaborators | -| Example 07 | Static inverse solve in a SPARC-like tokamak. | Anyone | -| Example 08 | Static inverse solve in an ITER-like tokamak. | Anyone | +| Example 07a | Static inverse solve in the "Anamak" toy tokamak. | Anyone | +| Example 07b | Static inverse solve in a SPARC-like tokamak. | Anyone | | Example 09 | Learn how to use and build virtual circuits for plasma shape control. | Anyone | | Example 10 | Learn how to calculate growth rates associated with vertically unstable modes. | Anyone | From c70e2244b8297833a59a1348671a6916e9c81729 Mon Sep 17 00:00:00 2001 From: kpentland Date: Thu, 26 Mar 2026 16:12:13 +0000 Subject: [PATCH 06/11] adding missing triangularity calculations --- examples/example07a - Anamak.ipynb | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/examples/example07a - Anamak.ipynb b/examples/example07a - Anamak.ipynb index 028e1535..f58fb73d 100644 --- a/examples/example07a - Anamak.ipynb +++ b/examples/example07a - Anamak.ipynb @@ -504,6 +504,16 @@ " )" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "print(f\"Triangularity = {eq.triangularity()}\")\n", + "print(f\"Avg. Z (Jtor) position = {eq.Zcurrent()} [m]\")" + ] + }, { "cell_type": "code", "execution_count": null, From 2af20ebf385c26c12f351b1230e8b4b1da1bd10c Mon Sep 17 00:00:00 2001 From: kpentland Date: Fri, 27 Mar 2026 12:03:29 +0000 Subject: [PATCH 07/11] fixing typo --- machine_configs/README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/machine_configs/README.md b/machine_configs/README.md index 5787cb39..8bcb3457 100644 --- a/machine_configs/README.md +++ b/machine_configs/README.md @@ -7,5 +7,5 @@ Here we store a number of **machine description** files in pickle format that de | test | A test tokamak used in the FreeGSNKE unit tests. | N\A | example | A simple toy tokamak. | Example0 notebook | MAST-U | A **MAST-U-like** tokamak. | UKAEA -| SPARC | A **SPARC-U-like** tokamak. | [Link](https://github.com/cfs-energy/SPARCPublic) +| SPARC | A **SPARC-like** tokamak. | [Link](https://github.com/cfs-energy/SPARCPublic) | ITER | An **ITER-like** tokamak. | [Link](https://github.com/ProjectTorreyPines/FUSE.jl) From 646b726883c4c14c41756fdd59435b018386755f Mon Sep 17 00:00:00 2001 From: kpentland Date: Mon, 13 Apr 2026 14:08:16 +0100 Subject: [PATCH 08/11] re-adding fixed ITER example from #67 --- examples/README.md | 1 + examples/example07c - ITER.ipynb | 319 +++++++++++++++++++++++++++++++ 2 files changed, 320 insertions(+) create mode 100644 examples/example07c - ITER.ipynb diff --git a/examples/README.md b/examples/README.md index ae1c2437..19f1cebe 100644 --- a/examples/README.md +++ b/examples/README.md @@ -14,6 +14,7 @@ These example Jupyter notebooks are intended to be the **first port of call for | Example 06a/b | Simulate (static) MAST-U equilibria over an entire shot using inputs from EFIT++ (requires internal UKAEA MAST-U database). | UKAEA employees + collaborators | | Example 07a | Static inverse solve in the "Anamak" toy tokamak. | Anyone | | Example 07b | Static inverse solve in a SPARC-like tokamak. | Anyone | +| Example 07c | Static inverse solve in an ITER-like tokamak. | Anyone | | Example 09 | Learn how to use and build virtual circuits for plasma shape control. | Anyone | | Example 10 | Learn how to calculate growth rates associated with vertically unstable modes. | Anyone | diff --git a/examples/example07c - ITER.ipynb b/examples/example07c - ITER.ipynb new file mode 100644 index 00000000..65091d8b --- /dev/null +++ b/examples/example07c - ITER.ipynb @@ -0,0 +1,319 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example: Static inverse free-boundary equilibrium calculations (in ITER)\n", + "\n", + "---\n", + "\n", + "Here we will generate an equilibrium (find coil currents with the inverse solver) in an ITER-like tokamak. \n", + "\n", + "The machine description comes from files located [here](https://github.com/ProjectTorreyPines/FUSE.jl).\n", + "\n", + "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create the machine object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build machine\n", + "from freegsnke import build_machine\n", + "tokamak = build_machine.tokamak(\n", + " active_coils_path=f\"../machine_configs/ITER/ITER_active_coils.pickle\",\n", + " passive_coils_path=f\"../machine_configs/ITER/ITER_passive_coils.pickle\",\n", + " limiter_path=f\"../machine_configs/ITER/ITER_limiter.pickle\",\n", + " wall_path=f\"../machine_configs/ITER/ITER_wall.pickle\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the machine\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "plt.tight_layout()\n", + "\n", + "tokamak.plot(axis=ax1, show=False)\n", + "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", + "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", + "\n", + "ax1.grid(alpha=0.5)\n", + "ax1.set_aspect('equal')\n", + "ax1.set_xlim(1.1, 12.4)\n", + "ax1.set_ylim(-8.3, 8.3)\n", + "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", + "ax1.set_ylabel(r'Height, $Z$ [m]')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate an equilibrium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import equilibrium_update\n", + "\n", + "eq = equilibrium_update.Equilibrium(\n", + " tokamak=tokamak, # provide tokamak object\n", + " Rmin=3.2, Rmax=8.8, # radial range\n", + " Zmin=-5, Zmax=5, # 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", + " # psi=plasma_psi\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate a profile object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke.jtor_update import Fiesta_Topeol\n", + "profiles = Fiesta_Topeol(\n", + " eq=eq,\n", + " Beta0=0.5978,\n", + " Ip=15e6,\n", + " fvac=0.5,\n", + " Raxis=6.2,\n", + " alpha_m=2.0,\n", + " alpha_n=1.395\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the static nonlinear solver" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import GSstaticsolver\n", + "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Constraints" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Rx = 5.02 # X-point radius\n", + "Zx = -3.23 # X-point height\n", + "\n", + "Ro = 6.34 # X-point radius\n", + "Zo = 0.66 # X-point height\n", + "\n", + "# set desired null_points locations\n", + "# this can include X-point and O-point locations\n", + "null_points = [[Rx, Ro], [Zx, Zo]]\n", + "\n", + "# set desired isoflux constraints with format \n", + "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", + "# with each isoflux_i = [R_coords, Z_coords]\n", + "isoflux_set = np.array([[\n", + " [4.25455147,\n", + " 4.1881875,\n", + " 4.2625,\n", + " 4.45683769,\n", + " 4.78746942,\n", + " 5.31835938,\n", + " 5.91875,\n", + " 6.44275166,\n", + " 6.92285156,\n", + " 7.35441013,\n", + " 7.73027344,\n", + " 8.03046875,\n", + " 8.19517497,\n", + " 8.05673414,\n", + " 7.75308013,\n", + " 7.37358451,\n", + " 6.94355469,\n", + " 6.47773438,\n", + " 5.99121094,\n", + " 5.49433594,\n", + " Rx,\n", + " 4.79042969,\n", + " 4.56269531,\n", + " 4.36038615], \n", + "[0.0,\n", + " 1.13554688,\n", + " 2.2565134,\n", + " 3.16757813,\n", + " 3.80507812,\n", + " 4.0499021,\n", + " 3.93089003,\n", + " 3.665625,\n", + " 3.31076604,\n", + " 2.86875,\n", + " 2.3199601,\n", + " 1.62832118,\n", + " 0.657421875,\n", + " -0.35859375,\n", + " -1.05585937,\n", + " -1.59375,\n", + " -2.03492727,\n", + " -2.41356403,\n", + " -2.75622016,\n", + " -3.08699278,\n", + " Zx,\n", + " -2.40548044,\n", + " -1.55982142,\n", + " -0.67734375]\n", + " ]])\n", + " \n", + "\n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(null_points=null_points,\n", + " isoflux_set=isoflux_set)\n", + "\n", + "\n", + "# remove from the coils available for control the radial field coil \n", + "# eq.tokamak['VS3'].control = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The inverse solve" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-4,\n", + " target_relative_psit_update=1e-3,\n", + " # max_solving_iterations=150,\n", + " # max_iter_per_update=10,\n", + " # max_rel_update_size=0.075,\n", + " # damping_factor=.99,\n", + " # max_rel_psit=.05,\n", + " verbose=True, # print output\n", + " l2_reg=1e-14,\n", + " Picard_handover=1e-4,\n", + " full_jacobian_handover=[1e-3, 1e-2]\n", + " )\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints\n", + "ax1.set_xlim(1.1, 12.4)\n", + "ax1.set_ylim(-8.3, 8.3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eq.tokamak.getCurrents()\n", + "\n", + "# # save coil currents to file\n", + "# import pickle\n", + "# with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", + "# pickle.dump(obj=inverse_current_values, file=f)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "freegsnke-dev", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.12" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From cd21223aa519b43c1fffadc5ccdd7b3c15936cc9 Mon Sep 17 00:00:00 2001 From: kpentland Date: Mon, 13 Apr 2026 14:10:36 +0100 Subject: [PATCH 09/11] adding fixed ITER example08, for removal later --- ...xample08 - static_inverse_solve_ITER.ipynb | 319 ++++++++++++++++++ 1 file changed, 319 insertions(+) create mode 100644 examples/example08 - static_inverse_solve_ITER.ipynb diff --git a/examples/example08 - static_inverse_solve_ITER.ipynb b/examples/example08 - static_inverse_solve_ITER.ipynb new file mode 100644 index 00000000..65091d8b --- /dev/null +++ b/examples/example08 - static_inverse_solve_ITER.ipynb @@ -0,0 +1,319 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Example: Static inverse free-boundary equilibrium calculations (in ITER)\n", + "\n", + "---\n", + "\n", + "Here we will generate an equilibrium (find coil currents with the inverse solver) in an ITER-like tokamak. \n", + "\n", + "The machine description comes from files located [here](https://github.com/ProjectTorreyPines/FUSE.jl).\n", + "\n", + "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Import packages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Create the machine object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# build machine\n", + "from freegsnke import build_machine\n", + "tokamak = build_machine.tokamak(\n", + " active_coils_path=f\"../machine_configs/ITER/ITER_active_coils.pickle\",\n", + " passive_coils_path=f\"../machine_configs/ITER/ITER_passive_coils.pickle\",\n", + " limiter_path=f\"../machine_configs/ITER/ITER_limiter.pickle\",\n", + " wall_path=f\"../machine_configs/ITER/ITER_wall.pickle\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the machine\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "plt.tight_layout()\n", + "\n", + "tokamak.plot(axis=ax1, show=False)\n", + "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", + "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", + "\n", + "ax1.grid(alpha=0.5)\n", + "ax1.set_aspect('equal')\n", + "ax1.set_xlim(1.1, 12.4)\n", + "ax1.set_ylim(-8.3, 8.3)\n", + "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", + "ax1.set_ylabel(r'Height, $Z$ [m]')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate an equilibrium" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import equilibrium_update\n", + "\n", + "eq = equilibrium_update.Equilibrium(\n", + " tokamak=tokamak, # provide tokamak object\n", + " Rmin=3.2, Rmax=8.8, # radial range\n", + " Zmin=-5, Zmax=5, # 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", + " # psi=plasma_psi\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Instantiate a profile object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke.jtor_update import Fiesta_Topeol\n", + "profiles = Fiesta_Topeol(\n", + " eq=eq,\n", + " Beta0=0.5978,\n", + " Ip=15e6,\n", + " fvac=0.5,\n", + " Raxis=6.2,\n", + " alpha_m=2.0,\n", + " alpha_n=1.395\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Load the static nonlinear solver" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from freegsnke import GSstaticsolver\n", + "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Constraints" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "Rx = 5.02 # X-point radius\n", + "Zx = -3.23 # X-point height\n", + "\n", + "Ro = 6.34 # X-point radius\n", + "Zo = 0.66 # X-point height\n", + "\n", + "# set desired null_points locations\n", + "# this can include X-point and O-point locations\n", + "null_points = [[Rx, Ro], [Zx, Zo]]\n", + "\n", + "# set desired isoflux constraints with format \n", + "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", + "# with each isoflux_i = [R_coords, Z_coords]\n", + "isoflux_set = np.array([[\n", + " [4.25455147,\n", + " 4.1881875,\n", + " 4.2625,\n", + " 4.45683769,\n", + " 4.78746942,\n", + " 5.31835938,\n", + " 5.91875,\n", + " 6.44275166,\n", + " 6.92285156,\n", + " 7.35441013,\n", + " 7.73027344,\n", + " 8.03046875,\n", + " 8.19517497,\n", + " 8.05673414,\n", + " 7.75308013,\n", + " 7.37358451,\n", + " 6.94355469,\n", + " 6.47773438,\n", + " 5.99121094,\n", + " 5.49433594,\n", + " Rx,\n", + " 4.79042969,\n", + " 4.56269531,\n", + " 4.36038615], \n", + "[0.0,\n", + " 1.13554688,\n", + " 2.2565134,\n", + " 3.16757813,\n", + " 3.80507812,\n", + " 4.0499021,\n", + " 3.93089003,\n", + " 3.665625,\n", + " 3.31076604,\n", + " 2.86875,\n", + " 2.3199601,\n", + " 1.62832118,\n", + " 0.657421875,\n", + " -0.35859375,\n", + " -1.05585937,\n", + " -1.59375,\n", + " -2.03492727,\n", + " -2.41356403,\n", + " -2.75622016,\n", + " -3.08699278,\n", + " Zx,\n", + " -2.40548044,\n", + " -1.55982142,\n", + " -0.67734375]\n", + " ]])\n", + " \n", + "\n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(null_points=null_points,\n", + " isoflux_set=isoflux_set)\n", + "\n", + "\n", + "# remove from the coils available for control the radial field coil \n", + "# eq.tokamak['VS3'].control = False" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### The inverse solve" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-4,\n", + " target_relative_psit_update=1e-3,\n", + " # max_solving_iterations=150,\n", + " # max_iter_per_update=10,\n", + " # max_rel_update_size=0.075,\n", + " # damping_factor=.99,\n", + " # max_rel_psit=.05,\n", + " verbose=True, # print output\n", + " l2_reg=1e-14,\n", + " Picard_handover=1e-4,\n", + " full_jacobian_handover=[1e-3, 1e-2]\n", + " )\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints\n", + "ax1.set_xlim(1.1, 12.4)\n", + "ax1.set_ylim(-8.3, 8.3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eq.tokamak.getCurrents()\n", + "\n", + "# # save coil currents to file\n", + "# import pickle\n", + "# with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", + "# pickle.dump(obj=inverse_current_values, file=f)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "freegsnke-dev", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.12" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} From 9e1c5abc806b0215fea40a21b911ccb912383444 Mon Sep 17 00:00:00 2001 From: kpentland Date: Mon, 13 Apr 2026 14:12:07 +0100 Subject: [PATCH 10/11] removing example08 after merging in 'main' --- ...xample08 - static_inverse_solve_ITER.ipynb | 319 ------------------ 1 file changed, 319 deletions(-) delete mode 100644 examples/example08 - static_inverse_solve_ITER.ipynb diff --git a/examples/example08 - static_inverse_solve_ITER.ipynb b/examples/example08 - static_inverse_solve_ITER.ipynb deleted file mode 100644 index 65091d8b..00000000 --- a/examples/example08 - static_inverse_solve_ITER.ipynb +++ /dev/null @@ -1,319 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# Example: Static inverse free-boundary equilibrium calculations (in ITER)\n", - "\n", - "---\n", - "\n", - "Here we will generate an equilibrium (find coil currents with the inverse solver) in an ITER-like tokamak. \n", - "\n", - "The machine description comes from files located [here](https://github.com/ProjectTorreyPines/FUSE.jl).\n", - "\n", - "The equilbirium\\profile parameters are **completely made up** - please experiment on your own and change them to more realistic values as you please!" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Import packages" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numpy as np" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Create the machine object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# build machine\n", - "from freegsnke import build_machine\n", - "tokamak = build_machine.tokamak(\n", - " active_coils_path=f\"../machine_configs/ITER/ITER_active_coils.pickle\",\n", - " passive_coils_path=f\"../machine_configs/ITER/ITER_passive_coils.pickle\",\n", - " limiter_path=f\"../machine_configs/ITER/ITER_limiter.pickle\",\n", - " wall_path=f\"../machine_configs/ITER/ITER_wall.pickle\",\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# plot the machine\n", - "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", - "plt.tight_layout()\n", - "\n", - "tokamak.plot(axis=ax1, show=False)\n", - "ax1.plot(tokamak.limiter.R, tokamak.limiter.Z, color='k', linewidth=1.2, linestyle=\"--\")\n", - "# ax1.plot(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, linestyle=\"-\")\n", - "\n", - "ax1.grid(alpha=0.5)\n", - "ax1.set_aspect('equal')\n", - "ax1.set_xlim(1.1, 12.4)\n", - "ax1.set_ylim(-8.3, 8.3)\n", - "ax1.set_xlabel(r'Major radius, $R$ [m]')\n", - "ax1.set_ylabel(r'Height, $Z$ [m]')" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Instantiate an equilibrium" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke import equilibrium_update\n", - "\n", - "eq = equilibrium_update.Equilibrium(\n", - " tokamak=tokamak, # provide tokamak object\n", - " Rmin=3.2, Rmax=8.8, # radial range\n", - " Zmin=-5, Zmax=5, # 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", - " # psi=plasma_psi\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Instantiate a profile object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke.jtor_update import Fiesta_Topeol\n", - "profiles = Fiesta_Topeol(\n", - " eq=eq,\n", - " Beta0=0.5978,\n", - " Ip=15e6,\n", - " fvac=0.5,\n", - " Raxis=6.2,\n", - " alpha_m=2.0,\n", - " alpha_n=1.395\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Load the static nonlinear solver" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from freegsnke import GSstaticsolver\n", - "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Constraints" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "Rx = 5.02 # X-point radius\n", - "Zx = -3.23 # X-point height\n", - "\n", - "Ro = 6.34 # X-point radius\n", - "Zo = 0.66 # X-point height\n", - "\n", - "# set desired null_points locations\n", - "# this can include X-point and O-point locations\n", - "null_points = [[Rx, Ro], [Zx, Zo]]\n", - "\n", - "# set desired isoflux constraints with format \n", - "# isoflux_set = [isoflux_0, isoflux_1 ... ] \n", - "# with each isoflux_i = [R_coords, Z_coords]\n", - "isoflux_set = np.array([[\n", - " [4.25455147,\n", - " 4.1881875,\n", - " 4.2625,\n", - " 4.45683769,\n", - " 4.78746942,\n", - " 5.31835938,\n", - " 5.91875,\n", - " 6.44275166,\n", - " 6.92285156,\n", - " 7.35441013,\n", - " 7.73027344,\n", - " 8.03046875,\n", - " 8.19517497,\n", - " 8.05673414,\n", - " 7.75308013,\n", - " 7.37358451,\n", - " 6.94355469,\n", - " 6.47773438,\n", - " 5.99121094,\n", - " 5.49433594,\n", - " Rx,\n", - " 4.79042969,\n", - " 4.56269531,\n", - " 4.36038615], \n", - "[0.0,\n", - " 1.13554688,\n", - " 2.2565134,\n", - " 3.16757813,\n", - " 3.80507812,\n", - " 4.0499021,\n", - " 3.93089003,\n", - " 3.665625,\n", - " 3.31076604,\n", - " 2.86875,\n", - " 2.3199601,\n", - " 1.62832118,\n", - " 0.657421875,\n", - " -0.35859375,\n", - " -1.05585937,\n", - " -1.59375,\n", - " -2.03492727,\n", - " -2.41356403,\n", - " -2.75622016,\n", - " -3.08699278,\n", - " Zx,\n", - " -2.40548044,\n", - " -1.55982142,\n", - " -0.67734375]\n", - " ]])\n", - " \n", - "\n", - "# instantiate the freegsnke constrain object\n", - "from freegsnke.inverse import Inverse_optimizer\n", - "constrain = Inverse_optimizer(null_points=null_points,\n", - " isoflux_set=isoflux_set)\n", - "\n", - "\n", - "# remove from the coils available for control the radial field coil \n", - "# eq.tokamak['VS3'].control = False" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### The inverse solve" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "GSStaticSolver.inverse_solve(eq=eq, \n", - " profiles=profiles, \n", - " constrain=constrain, \n", - " target_relative_tolerance=1e-4,\n", - " target_relative_psit_update=1e-3,\n", - " # max_solving_iterations=150,\n", - " # max_iter_per_update=10,\n", - " # max_rel_update_size=0.075,\n", - " # damping_factor=.99,\n", - " # max_rel_psit=.05,\n", - " verbose=True, # print output\n", - " l2_reg=1e-14,\n", - " Picard_handover=1e-4,\n", - " full_jacobian_handover=[1e-3, 1e-2]\n", - " )\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", - "\n", - "ax1.grid(zorder=0, alpha=0.75)\n", - "ax1.set_aspect('equal')\n", - "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", - "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", - "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", - "constrain.plot(axis=ax1, show=False) # plots the contraints\n", - "ax1.set_xlim(1.1, 12.4)\n", - "ax1.set_ylim(-8.3, 8.3)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "eq.tokamak.getCurrents()\n", - "\n", - "# # save coil currents to file\n", - "# import pickle\n", - "# with open('simple_diverted_currents_PaxisIp.pk', 'wb') as f:\n", - "# pickle.dump(obj=inverse_current_values, file=f)" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "freegsnke-dev", - "language": "python", - "name": "python3" - }, - "language_info": { - "codemirror_mode": { - "name": "ipython", - "version": 3 - }, - "file_extension": ".py", - "mimetype": "text/x-python", - "name": "python", - "nbconvert_exporter": "python", - "pygments_lexer": "ipython3", - "version": "3.13.12" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} From 2d7ddae7522f24608b0b6f2311174149d6725bb1 Mon Sep 17 00:00:00 2001 From: kpentland Date: Fri, 15 May 2026 13:52:38 +0100 Subject: [PATCH 11/11] Adding a snowflake example into the Anamak --- examples/example07a - Anamak.ipynb | 62 ++++++++++++++++++++++++++++-- 1 file changed, 58 insertions(+), 4 deletions(-) diff --git a/examples/example07a - Anamak.ipynb b/examples/example07a - Anamak.ipynb index f58fb73d..c7ed8486 100644 --- a/examples/example07a - Anamak.ipynb +++ b/examples/example07a - Anamak.ipynb @@ -750,19 +750,73 @@ "constrain.plot(axis=ax1, show=False) # plots the contraints" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now let us try to use the snowflake (i.e. second order null) constraints to make something a bit more exotic. " + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# null points constraints\n", + "Rx = 0.8 # X-point radius\n", + "Zx = 0.65 # X-point height\n", + "null_points_2nd_order = [[Rx, Rx], [Zx, -Zx]] \n", + "\n", + "# isoflux constraints (including null points here, not really needed but should be fine)\n", + "isoflux_set = np.array([\n", + " [\n", + " [Rx, Rx], \n", + " [Zx, -Zx]\n", + " ]]\n", + " )\n", + " \n", + "# instantiate the freegsnke constrain object\n", + "from freegsnke.inverse import Inverse_optimizer\n", + "constrain = Inverse_optimizer(\n", + " null_points_2nd_order=null_points_2nd_order,\n", + " # isoflux_set=isoflux_set,\n", + " )" + ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "# solve\n", + "GSStaticSolver.inverse_solve(eq=eq, \n", + " profiles=profiles, \n", + " constrain=constrain, \n", + " target_relative_tolerance=1e-5,\n", + " target_relative_psit_update=1e-3,\n", + " verbose=True,\n", + " l2_reg=np.array([1e-12]*N_actives),\n", + " force_up_down_symmetric=True,\n", + " )" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# plot the equilibrium\n", + "fig1, ax1 = plt.subplots(1, 1, figsize=(7, 15), dpi=80)\n", + "ax1.grid(zorder=0, alpha=0.75)\n", + "ax1.set_aspect('equal')\n", + "eq.tokamak.plot(axis=ax1,show=False) # plots the active coils and passive structures\n", + "ax1.fill(tokamak.wall.R, tokamak.wall.Z, color='k', linewidth=1.2, facecolor='w', zorder=0) # plots the limiter\n", + "eq.plot(axis=ax1,show=False) # plots the equilibrium\n", + "constrain.plot(axis=ax1, show=False) # plots the contraints" + ] }, { "cell_type": "code", @@ -774,7 +828,7 @@ ], "metadata": { "kernelspec": { - "display_name": "freegsnke4e", + "display_name": "freegsnke", "language": "python", "name": "python3" }, @@ -788,7 +842,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.16" + "version": "3.10.20" } }, "nbformat": 4,