Skip to content
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
58 changes: 36 additions & 22 deletions examples/example10 - growth_rates.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -246,6 +246,13 @@
"Here we will be using **all** available metal modes, however, later on we will show how to exclude those (via additional options to `nl_solver`) with specific frequencies, timescales, coupling, etc to reduce computational runtime (without losing too much information)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
Expand Down Expand Up @@ -350,7 +357,7 @@
"\n",
"# multiply each metal current (from the eigenvector) with its corresponding Greens matrix and sum\n",
"# (don't forget to omit the plasma current mode, i.e. the final element)\n",
"flux = np.sum(mode_currents[0:-1, np.newaxis, np.newaxis] * nonlinear_solver.vessel_modes_greens, axis=0)\n",
"flux = np.sum(mode_currents[0:-1, np.newaxis, np.newaxis] * nonlinear_solver.vessel_modes_greens[i], axis=0)\n",
"\n",
"# plot\n",
"fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n",
Expand All @@ -369,6 +376,13 @@
"cbar = plt.colorbar(im, ax=ax, fraction=0.09)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
Expand Down Expand Up @@ -503,7 +517,7 @@
"\n",
"This parameter quantifies the ability of the active coil system to stabilise vertical displacements of the plasma. It is purely inductive as it does not contain any resistance terms and is given by\n",
"\n",
"$$ m_s := \\lambda \\left[ -M_{m,m}^{-1} \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS} \\right) \\right], $$\n",
"$$ m_s := \\lambda \\left[ -M_{m,m}^{-1} \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} \\right) \\right], $$\n",
"\n",
"where the operator $\\lambda[A]$ returns the largest non-negative eigenvalue of the matrix $A$. This is the counterpart to the \"resistive\" stability margin (i.e. the deformable timescale) that is described in the \"Removing the plasma equation\" - where in this case the plasma equation has been omitted (i.e. $\\dot{I_p} = 0$).\n",
"\n",
Expand All @@ -530,7 +544,7 @@
"The (generalised) Leuer stability parameter is defined as the ratio of stabilising force gradient to de-stabilising force gradient acting on the system:\n",
"\n",
"\\begin{equation}\n",
" f := \\frac{f_{stab}}{f_{destab}} = \\frac{\\bold{I}_y^T M'_{y,m} M^{-1}_{m,m} M'_{m,y} \\bold{I}_y}{\\bold{I}_y^T M''_{y,m} \\bold{I}_m},\n",
" f := \\frac{f_{stab}}{f_{destab}} = \\frac{\\vec{I}_y^T M'_{y,m} M^{-1}_{m,m} M'_{m,y} \\vec{I}_y}{\\vec{I}_y^T M''_{y,m} \\vec{I}_m},\n",
"\\end{equation}\n",
"\n",
"where $M'$ and $M''$ denote the first and second derivatives of the mutual inductances with respect to the $Z$ coordinate, respectively. Other matrices are defined as they are above. We can find these derivatives using the relations:\n",
Expand Down Expand Up @@ -580,10 +594,10 @@
"source": [
"### Back to the deformable model: Removing the plasma equation\n",
"\n",
"The deformable growth rates can also be computed in the case where we assume $\\bold{\\dot{I}}_p = 0$. For this, the final row/column of the state space matrices above are omitted and the system becomes:\n",
"The deformable growth rates can also be computed in the case where we assume $\\vec{\\dot{I}}_p = 0$. For this, the final row/column of the state space matrices above are omitted and the system becomes:\n",
"\n",
"\\begin{equation*}\n",
" \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS} \\right) \\dot{\\bold{I}}_m + R_{m,m} \\bold{I}_m = \\bold{V}_m - M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{\\theta}}\\Bigg|_{GS} \\dot{\\bold{\\theta}},\n",
" \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} \\right) \\dot{\\vec{I}}_m + R_{m,m} \\vec{I}_m = \\vec{V}_m - M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}},\n",
"\\end{equation*}\n",
"\n",
"and we find the eigenvaluesas before. \n",
Expand Down Expand Up @@ -642,20 +656,20 @@
"\\end{equation*}\n",
"where $\\mathcal{I}_{act,act}$ is the identity matrix (with size equal to number of active coils). Supposing we retain $k$ columns of $P$ (including all of the active coils and a subset of the \"most strongly coupled\" passive structures), we can write the full vector of metal currents as \n",
"\\begin{equation*}\n",
" \\bold{I}_m = P \\bold{I}_d,\n",
" \\vec{I}_m = P \\vec{I}_d,\n",
"\\end{equation*}\n",
"where $\\bold{I}_d$ is the \"reduced\" set of $k$ mode currents.\n",
"where $\\vec{I}_d$ is the \"reduced\" set of $k$ mode currents.\n",
"\n",
"Now, recall the state space matrix system from above. If we apply our expression for $\\bold{I}_m$, we get:\n",
"Now, recall the state space matrix system from above. If we apply our expression for $\\vec{I}_m$, we get:\n",
"\n",
"\\begin{equation*}\n",
"\\begin{aligned}\n",
"\\left[\\begin{array}{cc}\n",
" M_{m,m}P + M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS}P & M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial I_p}\\Bigg|_{GS} \\\\ \n",
"\\frac{\\hat{\\bold{I}}_y}{R_p} \\left( M_{y,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{I}_m} \\Bigg|_{GS}P + M_{y,m}P \\right) & \\frac{\\hat{\\bold{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\bold{I}_y}{\\partial I_p}\\Bigg|_{GS}\n",
" M_{m,m}P + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_d} \\Bigg|_{GS} & M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\\\ \n",
"\\frac{\\hat{\\vec{I}}_y}{R_p} \\left( M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS}P + M_{y,m}P \\right) & \\frac{\\hat{\\vec{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS}\n",
"\\end{array}\\right]\n",
"\\left[\\begin{array}{c}\n",
"\\dot{\\bold{I}}_d \\\\ \n",
"\\dot{\\vec{I}}_d \\\\ \n",
"\\dot{I}_p\n",
"\\end{array}\\right] +\n",
"\n",
Expand All @@ -664,19 +678,19 @@
"0 & 1\n",
"\\end{array}\\right]\n",
"\\left[\\begin{array}{c}\n",
"\\bold{I}_d \\\\ \n",
"\\vec{I}_d \\\\ \n",
"I_p\n",
"\\end{array}\\right]\n",
"\n",
"&=\n",
"\\left[\\begin{array}{c}\n",
"\\bold{V}_m \\\\ \n",
"\\vec{V}_m \\\\ \n",
"0 \n",
"\\end{array}\\right]\n",
"-\n",
"\\left[\\begin{array}{c}\n",
"M_{m,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{\\theta}}\\Bigg|_{GS} \\dot{\\bold{\\theta}} \\\\ \n",
"\\frac{\\hat{\\bold{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\bold{I}_y}{\\partial \\bold{\\theta}}\\Bigg|_{GS} \\dot{\\bold{\\theta}}\n",
"M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} \\\\ \n",
"\\frac{\\hat{\\vec{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}}\n",
"\\end{array}\\right].\n",
"\n",
"\\end{aligned}\n",
Expand All @@ -698,9 +712,9 @@
"The process by which the passive modes are chosen is as follows:\n",
"1. A first \"estimate\" of\n",
" $$\n",
" \\lambda_i = \\left\\| \\left( \\frac{\\partial \\mathbf{I}_y}{\\partial \\mathbf{I}_m} \\right)_{:,i} \\right\\|,\n",
" \\lambda_i = \\left\\| \\left( \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\right)_{:,i} \\right\\|,\n",
" $$\n",
" is calculated for all passive structure modes (without fully solving the GS equation and calculating the true Jacobian). This is done by perturbing the currents in the modes (passive structures), updating the poloidal flux produced by these currents, and then evaluating the plasma current density $J_p(R,Z, \\psi)$ to obtain the perturbed $\\mathbf{I}_y$. Here, $\\psi$ is the total poloidal flux, i.e. the plasma and the metals contributions combined.\n",
" is calculated for all passive structure modes (without fully solving the GS equation and calculating the true Jacobian). This is done by perturbing the currents in the modes (passive structures), updating the poloidal flux produced by these currents, and then evaluating the plasma current density $J_p(R,Z, \\psi)$ to obtain the perturbed $\\vec{I}_y$. Here, $\\psi$ is the total poloidal flux, i.e. the plasma and the metals contributions combined.\n",
"\n",
"2. The $\\lambda_i$ are then sorted into ascending order and with $\\hat{\\lambda} = \\max_i \\{ \\lambda_i\\}$ being the largest (strongest coupling).\n",
"\n",
Expand Down Expand Up @@ -729,7 +743,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"Again, we can plot how these look using similar code as before, except we need to remember that now we've done a mode decomposition, the eigenmode currents $\\bold{I}_d$ need to be transformed back into the original metal currents $\\bold{I}_m$. This is done in the plotting cell below."
"Again, we can plot how these look using similar code as before, except we need to remember that now we've done a mode decomposition, the eigenmode currents $\\vec{I}_d$ need to be transformed back into the original metal currents $\\vec{I}_m$. This is done in the plotting cell below."
]
},
{
Expand Down Expand Up @@ -778,7 +792,7 @@
"print(f\"Timescale = {np.real(timescales[i]):.2e} [s]\")\n",
"\n",
"# multiply each metal current (from the eigenvector) with its corresponding Greens matrix and sum\n",
"flux = np.sum(mode_currents[:, np.newaxis, np.newaxis] * nonlinear_solver_option_1.vessel_modes_greens, axis=0)\n",
"flux = np.sum(mode_currents[:, np.newaxis, np.newaxis] * nonlinear_solver_option_1.vessel_modes_greens[i], axis=0)\n",
"\n",
"# plot\n",
"fig, ax = plt.subplots(1, 1, figsize=(5, 8), dpi=60)\n",
Expand Down Expand Up @@ -818,7 +832,7 @@
"\n",
"2. As in option 1, calculate a first estimate of \n",
" $$\n",
" \\lambda_i = \\left\\| \\left( \\frac{\\partial \\mathbf{I}_y}{\\partial \\mathbf{I}_m} \\right)_{:,i} \\right\\|\n",
" \\lambda_i = \\left\\| \\left( \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\right)_{:,i} \\right\\|\n",
" $$ \n",
" for all passive structure modes without fully solving the GS equation for each (this is done for flagged modes too). This is prior to the full Jacobian calculation later on.\n",
"\n",
Expand Down Expand Up @@ -1234,7 +1248,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "freegsnke",
"language": "python",
"name": "python3"
},
Expand All @@ -1248,7 +1262,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.16"
"version": "3.10.20"
}
},
"nbformat": 4,
Expand Down
Loading