{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Grid Basis Vectors\n",
"\n",
"## Introduction\n",
"\n",
"In this notebook, we visualize the basis vectors of the EMC3-EIRENE grid. The basis vectors are defined in the computational space, which is a curvilinear coordinate system.\n",
"There are two types of basis vectors: covariant and contravariant, each of which is defined by the following equation:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\\mathbf{b}_i &= \\frac{\\partial \\mathbf{x}}{\\partial q^i} \\\\\n",
"\\mathbf{b}^i &= \\frac{\\partial q^i}{\\partial \\mathbf{x}} = \\nabla q^i\n",
"\\end{align}\n",
"$$\n",
"\n",
"where $\\mathbf{x}$ is the position vector in the space, $q^i$ is the coordinate in the computational space, and $i$ is the index of the coordinate.\n",
"EMC3-EIRENE grids follow the curvilinear coordinates, and each $q^i$ denotes radial($\\rho$), poloidal($\\theta$), and toroidal($\\zeta$) coordinates, respectively.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from plotly import graph_objects as go\n",
"\n",
"from cherab.lhd.emc3 import CenterGrid, Grid"
]
},
{
"attachments": {},
"cell_type": "markdown",
"metadata": {},
"source": [
"## Plot center points of the grid\n",
"\n",
"In this section, we plot the center points of the grid in the computational space.\n",
"The center points are calculated with tretrahedral barycenters like:\n",
"\n",
"$$\n",
"\\mathrm{x}_{\\mathrm{center}} ≡ \\frac{1}{N}∑_{i=1}^{N}\\mathrm{x}_i\n",
"$$\n",
"\n",
"where $N$ is the number of tetrahedra involving one voxel (cell) and $\\mathrm{x}_i$ is the barycenter of the $i$-th tetrahedron.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# NOTE: Use the demo grid dataset which is openly available\n",
"dataset = \"grid-demo.nc\"\n",
"zone = \"zone0\"\n",
"index = \"coarse\"\n",
"\n",
"# Instantiate grid & center grid objects\n",
"grid = Grid(zone, dataset=dataset)\n",
"cg = CenterGrid(zone, index_type=index, dataset=dataset)\n",
"\n",
"# Get grid vertices in (x, y, z) coordinates\n",
"verts = grid.generate_vertices().reshape((*grid.shape, 3), order=\"F\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# create coarse's grid index\n",
"L, M, N = grid.shape\n",
"poloidal_step = 6\n",
"toroidal_step = 4\n",
"radial_indices = [i for i in range(0, L)]\n",
"poloidal_indices = [i for i in range(0, M, poloidal_step)]\n",
"toroidal_indices = [i for i in range(0, N, toroidal_step)]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def plot_grid(\n",
" fig: go.Figure,\n",
" range_radial: slice | None = None,\n",
" range_poloidal: slice | None = None,\n",
" range_toroidal: slice | None = None,\n",
"):\n",
" \"\"\"Plot the grid vertices in 3D.\n",
"\n",
" Parameters\n",
" ----------\n",
" fig : go.Figure\n",
" The figure to plot the grid vertices on.\n",
" range_radial : slice, optional\n",
" The radial range of the grid vertices to plot, by default is slice(0, 4).\n",
" range_poloidal : slice, optional\n",
" The poloidal range of the grid vertices to plot, by default is slice(0, 4).\n",
" range_toroidal : slice, optional\n",
" The toroidal range of the grid vertices to plot, by default is slice(0, 4).\n",
" \"\"\"\n",
" if range_radial is None:\n",
" range_radial = slice(1, 5)\n",
" if range_poloidal is None:\n",
" range_poloidal = slice(0, 4)\n",
" if range_toroidal is None:\n",
" range_toroidal = slice(0, 4)\n",
"\n",
" # plot grid vertices radially\n",
" for ζ_n in toroidal_indices[range_toroidal]:\n",
" for θ_m in poloidal_indices[range_poloidal]:\n",
" fig.add_trace(\n",
" go.Scatter3d(\n",
" x=verts[range_radial, θ_m, ζ_n, 0],\n",
" y=verts[range_radial, θ_m, ζ_n, 1],\n",
" z=verts[range_radial, θ_m, ζ_n, 2],\n",
" mode=\"lines\",\n",
" line=dict(color=\"blue\", width=1),\n",
" showlegend=False,\n",
" hoverinfo=\"none\",\n",
" )\n",
" )\n",
"\n",
" # plot grid vertices poloidally\n",
" for ζ_n in toroidal_indices[range_toroidal]:\n",
" for ρ_l in radial_indices[range_radial]:\n",
" fig.add_trace(\n",
" go.Scatter3d(\n",
" x=verts[ρ_l, poloidal_indices[range_poloidal], ζ_n, 0],\n",
" y=verts[ρ_l, poloidal_indices[range_poloidal], ζ_n, 1],\n",
" z=verts[ρ_l, poloidal_indices[range_poloidal], ζ_n, 2],\n",
" mode=\"lines\",\n",
" line=dict(color=\"blue\", width=1),\n",
" showlegend=False,\n",
" hoverinfo=\"none\",\n",
" )\n",
" )\n",
"\n",
" # plot grid vertices toroidally\n",
" for θ_m in poloidal_indices[range_poloidal]:\n",
" for ρ_l in radial_indices[range_radial]:\n",
" fig.add_trace(\n",
" go.Scatter3d(\n",
" x=verts[ρ_l, θ_m, toroidal_indices[range_toroidal], 0],\n",
" y=verts[ρ_l, θ_m, toroidal_indices[range_toroidal], 1],\n",
" z=verts[ρ_l, θ_m, toroidal_indices[range_toroidal], 2],\n",
" mode=\"lines\",\n",
" line=dict(color=\"blue\", width=1),\n",
" showlegend=False,\n",
" hoverinfo=\"none\",\n",
" )\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = go.Figure(\n",
" layout=go.Layout(\n",
" scene=dict(\n",
" xaxis=dict(visible=False),\n",
" yaxis=dict(visible=False),\n",
" zaxis=dict(visible=False),\n",
" # aspectmode=\"data\",\n",
" ),\n",
" margin=dict(l=0, r=0, b=0, t=0),\n",
" width=700,\n",
" height=500,\n",
" )\n",
")\n",
"\n",
"# plot grid range\n",
"range_radial_cg = slice(1, 4)\n",
"range_poloidal_cg = slice(0, 3)\n",
"range_toroidal_cg = slice(0, 3)\n",
"\n",
"# extract center grid coordinates within range and reshape to (N, 3)\n",
"cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"\n",
"# plot grid vertices\n",
"plot_grid(fig)\n",
"\n",
"# plot center points of cells\n",
"fig.add_trace(\n",
" go.Scatter3d(\n",
" x=cg_extracted[:, 0],\n",
" y=cg_extracted[:, 1],\n",
" z=cg_extracted[:, 2],\n",
" mode=\"markers\",\n",
" marker=dict(size=5, color=\"red\"),\n",
" showlegend=False,\n",
" hovertemplate=\"Center points
\" + \"x: %{x}
y: %{y}
z: %{z}
\"\n",
" \"\",\n",
" )\n",
")\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Covariant basis $\\mathbf{b}_i = \\frac{∂ \\mathbf{x}}{∂ q^i}$\n",
"\n",
"For example, the vector notation of $\\mathbf{b}_ρ$ can be expressed by $\\mathbf{b}_ρ = \\begin{bmatrix} x_ρ & y_ρ & z_ρ \\end{bmatrix}^\\mathsf{T}$, where $x_ρ ≡ ∂_ρ x$.\n",
"The partial difference is calculated with the finite difference method.\n",
"\n",
"Firstly, we calculate the distance between each cell along to each coordinate axis and generate coordinates values.\n",
"Afterwards we calculate the partial difference of the coordinates values with `numpy.gradient` function.\n",
"\n",
"These procedure is handled by `CurvCoords` class."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from cherab.lhd.emc3.curvilinear import CurvCoords\n",
"\n",
"grid = CurvCoords(cg)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def plot_bases(\n",
" fig,\n",
" bases,\n",
" centers,\n",
" sizeref=0.4,\n",
" colorscale=None,\n",
" name=\"Bases\",\n",
"):\n",
" \"\"\"Plot the bases of the grid in 3D.\n",
"\n",
" Parameters\n",
" ----------\n",
" fig : go.Figure\n",
" The figure to plot the bases on.\n",
" bases : (N, 3) ndarray\n",
" The bases of the grid in (x, y, z) coordinates.\n",
" centers : (N, 3) ndarray\n",
" The centers of the grid in (x, y, z) coordinates.\n",
" sizeref : float, optional\n",
" The size reference for the bases, by default 0.4.\n",
" colorscale : list, optional\n",
" The colorscale for the bases, by default [[0, \"rgb(255, 0, 0)\"], [1, \"rgb(255, 0, 0)\"]].\n",
" name : str, optional\n",
" The name of the bases, by default \"Bases\".\n",
" \"\"\"\n",
" fig.add_trace(\n",
" go.Cone(\n",
" x=centers[:, 0],\n",
" y=centers[:, 1],\n",
" z=centers[:, 2],\n",
" u=bases[:, 0],\n",
" v=bases[:, 1],\n",
" w=bases[:, 2],\n",
" sizemode=\"absolute\",\n",
" sizeref=sizeref,\n",
" anchor=\"tail\",\n",
" colorscale=colorscale,\n",
" showscale=False,\n",
" showlegend=True,\n",
" name=name,\n",
" text=name,\n",
" )\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Obtain covariant bases\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"b_rho = grid.b_rho\n",
"b_theta = grid.b_theta\n",
"b_zeta = grid.b_zeta"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Plot covariant basis $\\mathbf{b}_i$ with grids\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# extract covariant bases within the grid range and reshape to (N, 3)\n",
"br_extracted = b_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"bt_extracted = b_theta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"bz_extracted = b_zeta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"\n",
"# Create Figure object\n",
"fig = go.Figure(\n",
" layout=go.Layout(\n",
" scene=dict(\n",
" xaxis=dict(visible=False),\n",
" yaxis=dict(visible=False),\n",
" zaxis=dict(visible=False),\n",
" # aspectmode=\"data\",\n",
" ),\n",
" margin=dict(l=0, r=0, b=0, t=0),\n",
" width=700,\n",
" height=500,\n",
" legend=dict(yanchor=\"top\", y=0.99, xanchor=\"left\", x=0.01),\n",
" )\n",
")\n",
"\n",
"# plot grid\n",
"plot_grid(fig)\n",
"\n",
"# show covariant bases b_ρ\n",
"plot_bases(\n",
" fig,\n",
" br_extracted,\n",
" cg_extracted,\n",
" sizeref=0.4,\n",
" colorscale=[[0, \"rgb(255, 0, 0)\"], [1, \"rgb(255, 0, 0)\"]],\n",
" name=\"$\\\\mathbf{b}_ρ$\",\n",
")\n",
"\n",
"# show covariant bases b_θ\n",
"plot_bases(\n",
" fig,\n",
" bt_extracted,\n",
" cg_extracted,\n",
" sizeref=0.6,\n",
" colorscale=[[0, \"rgb(0, 255, 0)\"], [1, \"rgb(0, 255, 0)\"]],\n",
" name=\"$\\\\mathbf{b}_θ$\",\n",
")\n",
"\n",
"# show covariant bases b_ζ\n",
"plot_bases(\n",
" fig,\n",
" bz_extracted,\n",
" cg_extracted,\n",
" sizeref=0.2,\n",
" colorscale=[[0, \"rgb(0, 0, 255)\"], [1, \"rgb(0, 0, 255)\"]],\n",
" name=\"$\\\\mathbf{b}_ζ$\",\n",
")\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Contravariant basis $\\mathbf{b}^i = \\frac{∂ q^i}{∂ \\mathbf{x}} = ∇ q^i$\n",
"\n",
"The contravariant basis is calculated with covariant bases using the following equation:\n",
"\n",
"$$\n",
"\\begin{align}\n",
"\\mathbf{b}^i &= \\frac{\\mathbf{b}_j \\times \\mathbf{b}_k}{\\mathbf{b}_i \\cdot (\\mathbf{b}_j \\times \\mathbf{b}_k)}\n",
"\\end{align}\n",
"$$\n",
"\n",
"where $i, j, k$ are the cyclic permutation of $ρ, θ, ζ$.\n",
"\n",
"The donominator of the equation equals the determinant of the Jacobian matrix.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"b_sup_rho = grid.b_sup_rho\n",
"b_sup_theta = grid.b_sup_theta\n",
"b_sup_zeta = grid.b_sup_zeta"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# extract center grid coordinates and bases within range and reshape to (N, 3)\n",
"cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"br_extracted = b_sup_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"bt_extracted = b_sup_theta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape(\n",
" (-1, 3)\n",
")\n",
"bz_extracted = b_sup_zeta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"\n",
"# Create Figure object\n",
"fig = go.Figure(\n",
" layout=go.Layout(\n",
" scene=dict(\n",
" xaxis=dict(visible=False),\n",
" yaxis=dict(visible=False),\n",
" zaxis=dict(visible=False),\n",
" # aspectmode=\"data\",\n",
" ),\n",
" margin=dict(l=0, r=0, b=0, t=0),\n",
" width=700,\n",
" height=500,\n",
" legend=dict(yanchor=\"top\", y=0.99, xanchor=\"left\", x=0.01),\n",
" )\n",
")\n",
"\n",
"# plot grid\n",
"plot_grid(fig)\n",
"\n",
"# show contravariant bases b^ρ\n",
"plot_bases(\n",
" fig,\n",
" br_extracted,\n",
" cg_extracted,\n",
" sizeref=0.3,\n",
" colorscale=[[0, \"rgb(255, 0, 255)\"], [1, \"rgb(255, 0, 255)\"]],\n",
" name=\"$\\\\mathbf{b}^ρ$\",\n",
")\n",
"\n",
"# show contravariant bases b^θ\n",
"plot_bases(\n",
" fig,\n",
" bt_extracted,\n",
" cg_extracted,\n",
" sizeref=0.5,\n",
" colorscale=[[0, \"rgb(255, 155, 0)\"], [1, \"rgb(255, 155, 0)\"]],\n",
" name=\"$\\\\mathbf{b}^θ$\",\n",
")\n",
"\n",
"# show contravariant bases b^ζ\n",
"plot_bases(\n",
" fig,\n",
" bz_extracted,\n",
" cg_extracted,\n",
" sizeref=0.2,\n",
" colorscale=[[0, \"rgb(0, 255, 255)\"], [1, \"rgb(0, 255, 255)\"]],\n",
" name=\"$\\\\mathbf{b}^ζ$\",\n",
")\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Compare co/contra-variant bases\n",
"\n",
"In this section, we compare the covariant and contravariant bases.\n",
"They are orthogonal to each other, and the inner product of them is equal to 1, i.e.\n",
"\n",
"$$\n",
"\\mathbf{b}^i \\cdot \\mathbf{b}_j = \\mathbf{b}_i \\cdot \\mathbf{b}^j = \\delta^i_j,\n",
"$$\n",
"\n",
"where $\\delta^i_j$ is the Kronecker delta."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $\\mathbf{b}^θ$ vs $\\mathbf{b}_ρ$"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = go.Figure(\n",
" layout=go.Layout(\n",
" scene=dict(\n",
" xaxis=dict(visible=False),\n",
" yaxis=dict(visible=False),\n",
" zaxis=dict(visible=False),\n",
" aspectmode=\"data\",\n",
" ),\n",
" margin=dict(l=0, r=0, b=0, t=0),\n",
" width=700,\n",
" height=500,\n",
" legend=dict(yanchor=\"top\", y=0.99, xanchor=\"left\", x=0.01),\n",
" )\n",
")\n",
"\n",
"# plot grid range\n",
"range_radial, range_radial_cg = slice(8, 12), slice(8, 11)\n",
"range_poloidal, range_poloidal_cg = slice(0, 2), slice(0, 1)\n",
"range_toroidal, range_toroidal_cg = slice(0, 3), slice(0, 2)\n",
"\n",
"# extract center grid coordinates and covariant bases within range and reshape to (N, 3)\n",
"cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"br_extracted = b_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"b_sup_t_extracted = b_sup_theta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape(\n",
" (-1, 3)\n",
")\n",
"\n",
"# plot grid\n",
"plot_grid(fig, range_radial, range_poloidal, range_toroidal)\n",
"\n",
"# show covariant bases b_ρ\n",
"plot_bases(\n",
" fig,\n",
" br_extracted,\n",
" cg_extracted,\n",
" sizeref=0.45,\n",
" colorscale=[[0, \"rgb(255, 0, 0)\"], [1, \"rgb(255, 0, 0)\"]],\n",
" name=\"$\\\\mathbf{b}_ρ$\",\n",
")\n",
"\n",
"# show contravariant bases b^θ\n",
"plot_bases(\n",
" fig,\n",
" b_sup_t_extracted,\n",
" cg_extracted,\n",
" sizeref=0.6,\n",
" colorscale=[[0, \"rgb(255, 155, 0)\"], [1, \"rgb(255, 155, 0)\"]],\n",
" name=\"$\\\\mathbf{b}^θ$\",\n",
")\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### $\\mathbf{b}^\\rho$ vs $\\mathbf{b}_\\rho$\n",
"\n",
"In not orthogonal coordinates, $\\mathbf{b}^i$ is not equal to $\\mathbf{b}_i$."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"fig = go.Figure(\n",
" layout=go.Layout(\n",
" scene=dict(\n",
" xaxis=dict(visible=False),\n",
" yaxis=dict(visible=False),\n",
" zaxis=dict(visible=False),\n",
" aspectmode=\"data\",\n",
" ),\n",
" margin=dict(l=0, r=0, b=0, t=0),\n",
" width=700,\n",
" height=500,\n",
" legend=dict(yanchor=\"top\", y=0.99, xanchor=\"left\", x=0.01),\n",
" )\n",
")\n",
"\n",
"# plot grid range\n",
"range_radial, range_radial_cg = slice(8, 12), slice(8, 11)\n",
"range_poloidal, range_poloidal_cg = slice(0, 2), slice(0, 1)\n",
"range_toroidal, range_toroidal_cg = slice(0, 3), slice(0, 2)\n",
"\n",
"# extract center grid coordinates and covariant bases within range and reshape to (N, 3)\n",
"cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"br_extracted = b_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))\n",
"b_sup_r_extracted = b_sup_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape(\n",
" (-1, 3)\n",
")\n",
"\n",
"# plot grid\n",
"plot_grid(fig, range_radial, range_poloidal, range_toroidal)\n",
"\n",
"# show covariant bases b_ρ\n",
"plot_bases(\n",
" fig,\n",
" br_extracted,\n",
" cg_extracted,\n",
" sizeref=0.4,\n",
" colorscale=[[0, \"rgb(255, 0, 0)\"], [1, \"rgb(255, 0, 0)\"]],\n",
" name=\"$\\\\mathbf{b}_ρ$\",\n",
")\n",
"\n",
"# show contravariant bases b^ρ\n",
"plot_bases(\n",
" fig,\n",
" b_sup_r_extracted,\n",
" cg_extracted,\n",
" sizeref=0.4,\n",
" colorscale=[[0, \"rgb(255, 0, 255)\"], [1, \"rgb(255, 0, 255)\"]],\n",
" name=\"$\\\\mathbf{b}^ρ$\",\n",
")\n",
"\n",
"fig.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"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.12.10"
}
},
"nbformat": 4,
"nbformat_minor": 4
}