Note

This page was generated from notebooks/emc3/curvilinear_bases.ipynb.

Grid Basis Vectors

Introduction

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. There are two types of basis vectors: covariant and contravariant, each of which is defined by the following equation:

\[\begin{split}\begin{align} \mathbf{b}_i &= \frac{\partial \mathbf{x}}{\partial q^i} \\ \mathbf{b}^i &= \frac{\partial q^i}{\partial \mathbf{x}} = \nabla q^i \end{align}\end{split}\]

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. EMC3-EIRENE grids follow the curvilinear coordinates, and each \(q^i\) denotes radial(\(\rho\)), poloidal(\(\theta\)), and toroidal(\(\zeta\)) coordinates, respectively.

from plotly import graph_objects as go

from cherab.lhd.emc3 import CenterGrid, Grid

Plot center points of the grid

In this section, we plot the center points of the grid in the computational space. The center points are calculated with tretrahedral barycenters like:

\[\mathrm{x}_{\mathrm{center}} ≡ \frac{1}{N}∑_{i=1}^{N}\mathrm{x}_i\]

where \(N\) is the number of tetrahedra involving one voxel (cell) and \(\mathrm{x}_i\) is the barycenter of the \(i\)-th tetrahedron.

# NOTE: Use the demo grid dataset which is openly available
dataset = "grid-demo.nc"
zone = "zone0"
index = "coarse"

# Instantiate grid & center grid objects
grid = Grid(zone, dataset=dataset)
cg = CenterGrid(zone, index_type=index, dataset=dataset)

# Get grid vertices in (x, y, z) coordinates
verts = grid.generate_vertices().reshape((*grid.shape, 3), order="F")

Downloading file 'grid-demo.nc' from 'doi:10.5281/zenodo.15412887/grid-demo.nc' to '/home/docs/.cache/cherab/lhd'.

# create coarse's grid index
L, M, N = grid.shape
poloidal_step = 6
toroidal_step = 4
radial_indices = [i for i in range(0, L)]
poloidal_indices = [i for i in range(0, M, poloidal_step)]
toroidal_indices = [i for i in range(0, N, toroidal_step)]

def plot_grid(
    fig: go.Figure,
    range_radial: slice | None = None,
    range_poloidal: slice | None = None,
    range_toroidal: slice | None = None,
):
    """Plot the grid vertices in 3D.

    Parameters
    ----------
    fig : go.Figure
        The figure to plot the grid vertices on.
    range_radial : slice, optional
        The radial range of the grid vertices to plot, by default is slice(0, 4).
    range_poloidal : slice, optional
        The poloidal range of the grid vertices to plot, by default is slice(0, 4).
    range_toroidal : slice, optional
        The toroidal range of the grid vertices to plot, by default is slice(0, 4).
    """
    if range_radial is None:
        range_radial = slice(1, 5)
    if range_poloidal is None:
        range_poloidal = slice(0, 4)
    if range_toroidal is None:
        range_toroidal = slice(0, 4)

    # plot grid vertices radially
    for ζ_n in toroidal_indices[range_toroidal]:
        for θ_m in poloidal_indices[range_poloidal]:
            fig.add_trace(
                go.Scatter3d(
                    x=verts[range_radial, θ_m, ζ_n, 0],
                    y=verts[range_radial, θ_m, ζ_n, 1],
                    z=verts[range_radial, θ_m, ζ_n, 2],
                    mode="lines",
                    line=dict(color="blue", width=1),
                    showlegend=False,
                    hoverinfo="none",
                )
            )

    # plot grid vertices poloidally
    for ζ_n in toroidal_indices[range_toroidal]:
        for ρ_l in radial_indices[range_radial]:
            fig.add_trace(
                go.Scatter3d(
                    x=verts[ρ_l, poloidal_indices[range_poloidal], ζ_n, 0],
                    y=verts[ρ_l, poloidal_indices[range_poloidal], ζ_n, 1],
                    z=verts[ρ_l, poloidal_indices[range_poloidal], ζ_n, 2],
                    mode="lines",
                    line=dict(color="blue", width=1),
                    showlegend=False,
                    hoverinfo="none",
                )
            )

    # plot grid vertices toroidally
    for θ_m in poloidal_indices[range_poloidal]:
        for ρ_l in radial_indices[range_radial]:
            fig.add_trace(
                go.Scatter3d(
                    x=verts[ρ_l, θ_m, toroidal_indices[range_toroidal], 0],
                    y=verts[ρ_l, θ_m, toroidal_indices[range_toroidal], 1],
                    z=verts[ρ_l, θ_m, toroidal_indices[range_toroidal], 2],
                    mode="lines",
                    line=dict(color="blue", width=1),
                    showlegend=False,
                    hoverinfo="none",
                )
            )

fig = go.Figure(
    layout=go.Layout(
        scene=dict(
            xaxis=dict(visible=False),
            yaxis=dict(visible=False),
            zaxis=dict(visible=False),
            # aspectmode="data",
        ),
        margin=dict(l=0, r=0, b=0, t=0),
        width=700,
        height=500,
    )
)

# plot grid range
range_radial_cg = slice(1, 4)
range_poloidal_cg = slice(0, 3)
range_toroidal_cg = slice(0, 3)

# extract center grid coordinates within range and reshape to (N, 3)
cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))

# plot grid vertices
plot_grid(fig)

# plot center points of cells
fig.add_trace(
    go.Scatter3d(
        x=cg_extracted[:, 0],
        y=cg_extracted[:, 1],
        z=cg_extracted[:, 2],
        mode="markers",
        marker=dict(size=5, color="red"),
        showlegend=False,
        hovertemplate="<b>Center points</b><br>" + "x: %{x}<br>y: %{y}<br>z: %{z}<br>"
        "<extra></extra>",
    )
)

fig.show()

Covariant basis \(\mathbf{b}_i = \frac{∂ \mathbf{x}}{∂ q^i}\)

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\). The partial difference is calculated with the finite difference method.

Firstly, we calculate the distance between each cell along to each coordinate axis and generate coordinates values. Afterwards we calculate the partial difference of the coordinates values with numpy.gradient function.

These procedure is handled by CurvCoords class.

from cherab.lhd.emc3.curvilinear import CurvCoords

grid = CurvCoords(cg)

def plot_bases(
    fig,
    bases,
    centers,
    sizeref=0.4,
    colorscale=None,
    name="Bases",
):
    """Plot the bases of the grid in 3D.

    Parameters
    ----------
    fig : go.Figure
        The figure to plot the bases on.
    bases : (N, 3) ndarray
        The bases of the grid in (x, y, z) coordinates.
    centers : (N, 3) ndarray
        The centers of the grid in (x, y, z) coordinates.
    sizeref : float, optional
        The size reference for the bases, by default 0.4.
    colorscale : list, optional
        The colorscale for the bases, by default [[0, "rgb(255, 0, 0)"], [1, "rgb(255, 0, 0)"]].
    name : str, optional
        The name of the bases, by default "Bases".
    """
    fig.add_trace(
        go.Cone(
            x=centers[:, 0],
            y=centers[:, 1],
            z=centers[:, 2],
            u=bases[:, 0],
            v=bases[:, 1],
            w=bases[:, 2],
            sizemode="absolute",
            sizeref=sizeref,
            anchor="tail",
            colorscale=colorscale,
            showscale=False,
            showlegend=True,
            name=name,
            text=name,
        )
    )

Obtain covariant bases

b_rho = grid.b_rho
b_theta = grid.b_theta
b_zeta = grid.b_zeta

Plot covariant basis \(\mathbf{b}_i\) with grids

# extract covariant bases within the grid range and reshape to (N, 3)
br_extracted = b_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
bt_extracted = b_theta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
bz_extracted = b_zeta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))

# Create Figure object
fig = go.Figure(
    layout=go.Layout(
        scene=dict(
            xaxis=dict(visible=False),
            yaxis=dict(visible=False),
            zaxis=dict(visible=False),
            # aspectmode="data",
        ),
        margin=dict(l=0, r=0, b=0, t=0),
        width=700,
        height=500,
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
    )
)

# plot grid
plot_grid(fig)

# show covariant bases b_ρ
plot_bases(
    fig,
    br_extracted,
    cg_extracted,
    sizeref=0.4,
    colorscale=[[0, "rgb(255, 0, 0)"], [1, "rgb(255, 0, 0)"]],
    name="$\\mathbf{b}_ρ$",
)

# show covariant bases b_θ
plot_bases(
    fig,
    bt_extracted,
    cg_extracted,
    sizeref=0.6,
    colorscale=[[0, "rgb(0, 255, 0)"], [1, "rgb(0, 255, 0)"]],
    name="$\\mathbf{b}_θ$",
)

# show covariant bases b_ζ
plot_bases(
    fig,
    bz_extracted,
    cg_extracted,
    sizeref=0.2,
    colorscale=[[0, "rgb(0, 0, 255)"], [1, "rgb(0, 0, 255)"]],
    name="$\\mathbf{b}_ζ$",
)

fig.show()

Contravariant basis \(\mathbf{b}^i = \frac{∂ q^i}{∂ \mathbf{x}} = ∇ q^i\)

The contravariant basis is calculated with covariant bases using the following equation:

\[\begin{align} \mathbf{b}^i &= \frac{\mathbf{b}_j \times \mathbf{b}_k}{\mathbf{b}_i \cdot (\mathbf{b}_j \times \mathbf{b}_k)} \end{align}\]

where \(i, j, k\) are the cyclic permutation of \(ρ, θ, ζ\).

The donominator of the equation equals the determinant of the Jacobian matrix.

b_sup_rho = grid.b_sup_rho
b_sup_theta = grid.b_sup_theta
b_sup_zeta = grid.b_sup_zeta

# extract center grid coordinates and bases within range and reshape to (N, 3)
cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
br_extracted = b_sup_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
bt_extracted = b_sup_theta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape(
    (-1, 3)
)
bz_extracted = b_sup_zeta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))

# Create Figure object
fig = go.Figure(
    layout=go.Layout(
        scene=dict(
            xaxis=dict(visible=False),
            yaxis=dict(visible=False),
            zaxis=dict(visible=False),
            # aspectmode="data",
        ),
        margin=dict(l=0, r=0, b=0, t=0),
        width=700,
        height=500,
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
    )
)

# plot grid
plot_grid(fig)

# show contravariant bases b^ρ
plot_bases(
    fig,
    br_extracted,
    cg_extracted,
    sizeref=0.3,
    colorscale=[[0, "rgb(255, 0, 255)"], [1, "rgb(255, 0, 255)"]],
    name="$\\mathbf{b}^ρ$",
)

# show contravariant bases b^θ
plot_bases(
    fig,
    bt_extracted,
    cg_extracted,
    sizeref=0.5,
    colorscale=[[0, "rgb(255, 155, 0)"], [1, "rgb(255, 155, 0)"]],
    name="$\\mathbf{b}^θ$",
)

# show contravariant bases b^ζ
plot_bases(
    fig,
    bz_extracted,
    cg_extracted,
    sizeref=0.2,
    colorscale=[[0, "rgb(0, 255, 255)"], [1, "rgb(0, 255, 255)"]],
    name="$\\mathbf{b}^ζ$",
)

fig.show()

Compare co/contra-variant bases

In this section, we compare the covariant and contravariant bases. They are orthogonal to each other, and the inner product of them is equal to 1, i.e.

\[\mathbf{b}^i \cdot \mathbf{b}_j = \mathbf{b}_i \cdot \mathbf{b}^j = \delta^i_j,\]

where \(\delta^i_j\) is the Kronecker delta.

\(\mathbf{b}^θ\) vs \(\mathbf{b}_ρ\)

fig = go.Figure(
    layout=go.Layout(
        scene=dict(
            xaxis=dict(visible=False),
            yaxis=dict(visible=False),
            zaxis=dict(visible=False),
            aspectmode="data",
        ),
        margin=dict(l=0, r=0, b=0, t=0),
        width=700,
        height=500,
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
    )
)

# plot grid range
range_radial, range_radial_cg = slice(8, 12), slice(8, 11)
range_poloidal, range_poloidal_cg = slice(0, 2), slice(0, 1)
range_toroidal, range_toroidal_cg = slice(0, 3), slice(0, 2)

# extract center grid coordinates and covariant bases within range and reshape to (N, 3)
cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
br_extracted = b_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
b_sup_t_extracted = b_sup_theta[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape(
    (-1, 3)
)

# plot grid
plot_grid(fig, range_radial, range_poloidal, range_toroidal)

# show covariant bases b_ρ
plot_bases(
    fig,
    br_extracted,
    cg_extracted,
    sizeref=0.45,
    colorscale=[[0, "rgb(255, 0, 0)"], [1, "rgb(255, 0, 0)"]],
    name="$\\mathbf{b}_ρ$",
)

# show contravariant bases b^θ
plot_bases(
    fig,
    b_sup_t_extracted,
    cg_extracted,
    sizeref=0.6,
    colorscale=[[0, "rgb(255, 155, 0)"], [1, "rgb(255, 155, 0)"]],
    name="$\\mathbf{b}^θ$",
)

fig.show()

\(\mathbf{b}^\rho\) vs \(\mathbf{b}_\rho\)

In not orthogonal coordinates, \(\mathbf{b}^i\) is not equal to \(\mathbf{b}_i\).

fig = go.Figure(
    layout=go.Layout(
        scene=dict(
            xaxis=dict(visible=False),
            yaxis=dict(visible=False),
            zaxis=dict(visible=False),
            aspectmode="data",
        ),
        margin=dict(l=0, r=0, b=0, t=0),
        width=700,
        height=500,
        legend=dict(yanchor="top", y=0.99, xanchor="left", x=0.01),
    )
)

# plot grid range
range_radial, range_radial_cg = slice(8, 12), slice(8, 11)
range_poloidal, range_poloidal_cg = slice(0, 2), slice(0, 1)
range_toroidal, range_toroidal_cg = slice(0, 3), slice(0, 2)

# extract center grid coordinates and covariant bases within range and reshape to (N, 3)
cg_extracted = cg[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
br_extracted = b_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape((-1, 3))
b_sup_r_extracted = b_sup_rho[range_radial_cg, range_poloidal_cg, range_toroidal_cg, :].reshape(
    (-1, 3)
)

# plot grid
plot_grid(fig, range_radial, range_poloidal, range_toroidal)

# show covariant bases b_ρ
plot_bases(
    fig,
    br_extracted,
    cg_extracted,
    sizeref=0.4,
    colorscale=[[0, "rgb(255, 0, 0)"], [1, "rgb(255, 0, 0)"]],
    name="$\\mathbf{b}_ρ$",
)

# show contravariant bases b^ρ
plot_bases(
    fig,
    b_sup_r_extracted,
    cg_extracted,
    sizeref=0.4,
    colorscale=[[0, "rgb(255, 0, 255)"], [1, "rgb(255, 0, 255)"]],
    name="$\\mathbf{b}^ρ$",
)

fig.show()