Flux Surfaces
try:
import SBGeom
import SBGeom.VMEC as VMEC
import SBGeom.Mesh as Mesh
except:
import StellBlanket.SBGeom as SBGeom
import StellBlanket.SBGeom.VMEC as VMEC
import StellBlanket.SBGeom.Mesh as Mesh
import numpy as np
import plotly.graph_objects as go
import plotly
import traceback
plotly.offline.init_notebook_mode()
from plotly.subplots import make_subplots
Flux Surfaces Basics¶
The VMEC representation uses arrays $R^s_{mn}, Z^s_{mn}$ for each surface $s$ and calculates the position as a function of poloidal angle $\theta$ and toroidal angle $\phi$ as:
$$ R(s, \theta, \phi) = \sum_{mn} R^s_{mn} \cos(m \theta - n S \phi) $$ $$ Z(s, \theta, \phi) = \sum_{mn} R^s_{mn} \sin(m \theta - n S \phi) $$
where $S$ is the symmetry of the stellarator.
In order to parametrically define a stellarator, it is useful to have a notion of distance beyond the Last-Closed Flux Surface (LCFS). Here, this is done by considering the normal distance to the LCFS.
We can compute the (unnormalized) normal vector on the last closed flux surface in cylindrical coordinates $(R,Z,\phi)$:
$$ \vec{n}(s,\theta,\phi) = \frac{d\vec{r}}{d\theta} \times \frac{d\vec{r}}{d\phi} = \begin{bmatrix}\frac{\partial R(1, \theta,\phi)}{\partial\theta} \\ \frac{\partial Z(1, \theta,\phi)}{\partial\theta} \\ 0 \end{bmatrix} \times \begin{bmatrix}\frac{\partial R(1, \theta,\phi)}{\partial\phi} \\ \frac{\partial Z(1, \theta,\phi)}{\partial\phi} \\ -R \end{bmatrix} = \begin{bmatrix} - R\frac{ \partial Z}{\partial\theta} \\ R\frac{ \partial R}{\partial\theta} \\ \frac{\partial R}{\partial \theta} \frac{\partial Z }{\partial \phi} - \frac{\partial R}{\partial \phi} \frac{\partial Z }{\partial \theta} \end{bmatrix} $$ where $-R$ arises because $\partial\vec{e}_\phi/ \partial \phi = -R$. Computing this normal vector, normalizing it, multiplying it by a value $d$ and adding it to the point $(R(\theta,\phi),Z(\theta, \phi),\phi)$ yields a point a distance $d$ from the last closed flux surface.
In the code, we can load a VMEC file from an HDF5 as follows:
fs = VMEC.Flux_Surfaces_From_HDF5("helias5_vmec.nc4")
SYMM = fs.flux_surface_settings().symmetry
However, this does not have the capability of extending beyond the LCFS yet:
try:
fs.Return_Position(s=1.0, LCFS_distance_label = 1.0, theta = 0.5 * np.pi, phi = 0.0)
except Exception as e:
traceback.print_exc()
Traceback (most recent call last):
File "/tmp/ipykernel_52961/2877269816.py", line 2, in <module>
fs.Return_Position(s=1.0, LCFS_distance_label = 1.0, theta = 0.5 * np.pi, phi = 0.0)
~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
ValueError: Trying to have a distance from the LCFS with the base Flux_Surfaces class.
Therefore, we create a Flux_Surfaces_Normal_Extended which explicitly uses the normal vector to extend:
fs_norm = SBGeom.Flux_Surfaces_Normal_Extended(fs)
try:
print("Succes! Result: ", fs_norm.Return_Position(s=1.0, LCFS_distance_label=1.0, theta= 0.5 * np.pi, phi=0.5 * np.pi))
except Exception as e:
traceback.print_exc()
Succes! Result: [-0.16024207 21.90526585 4.29113531]
We can use this new object to create a surface a specific distance from the LCFS:
s = 1.0
d = 1.4
phi_min = 0.0
phi_max = 2 * np.pi / ( 2 * SYMM ) # Half Module is half of the symmetry due to Stellarator Symmetry.
ntheta = 50
nphi = 60
surface = fs_norm.Mesh_Surface(s = s, LCFS_distance_label = d, N_lines_theta = ntheta, N_lines_phi = nphi, phi_min = phi_min, phi_max = phi_max)
lcfs = fs_norm.Mesh_Surface(s = s, LCFS_distance_label = 0.0, N_lines_theta = ntheta, N_lines_phi = nphi, phi_min = phi_min - 0.5 * (phi_max - phi_min), phi_max = phi_max + 0.5 * (phi_max - phi_min))
To visualize, we use the Plotly package, which allows for an interactive, three-dimensional plot. We will also plot the planes on the edges of the desired $\phi$.
#====================================================================================================================================================================================================
# Plotting:
#====================================================================================================================================================================================================
def Plot_Phi_Plane(phi, fig, rmin,rmax, zmin,zmax, alpha = 1, color = 'blue', **kwargs):
'''
Plots a simple plane
'''
x = [rmin * np.cos(phi), rmax * np.cos(phi), rmax * np.cos(phi), rmin * np.cos(phi), rmin * np.cos(phi)]
y = [rmin * np.sin(phi), rmax * np.sin(phi), rmax * np.sin(phi), rmin * np.sin(phi), rmin * np.sin(phi)]
z = [zmin, zmin, zmax, zmax, zmin]
fig = fig.add_trace(go.Scatter3d(x=x, y=y, z=z, surfaceaxis=1, opacity=alpha, surfacecolor=color, mode='lines', line=dict(color='black', width=1), showlegend=False, name="phi plane"), **kwargs)
fig = fig.add_trace(go.Scatter3d(x=x, y=y, z=z, mode='lines', line=dict(color='black', width=3), showlegend=False), **kwargs)
def create_figure_planes():
fig = go.Figure()
fig.update_layout(autosize=False, width=600, height=500, margin=dict(l=0, r=0, b=0, t=0))
Plot_Phi_Plane(phi = phi_min, fig = fig, rmin = 16, rmax=28, zmin=-6, zmax=6, alpha = 0.7, color = 'lightblue')
Plot_Phi_Plane(phi = phi_max, fig = fig, rmin = 16, rmax=28, zmin=-6, zmax=6, alpha = 0.7, color = 'lightblue')
return fig
fig = create_figure_planes()
Mesh.Plot(fig, surface, color = 'red', name= "Normal Vector")
Mesh.Plot(fig, lcfs, color = 'rgb(200, 130, 200)', name= "LCFS", opacity = 0.8)
fig.show()
However, as can be seen, although we only evaluate from $\phi_{\text{min}}$ to $\phi_{\text{max}}$, the surface crosses the $\phi$ planes! This is because the equation of the normal vector above has a finite $\phi$ component. This crossing is undesirable for two reasons:
- Periodic (cylindrical) boundary conditions cannot be applied to the edges of a full module
- When transforming to a Fourier parametrisation of the blanket, the points are no longer equidistantly and uniformly spaced in $\phi$
To solve these issues, a constant $\phi$ version is also available. This performs several Newton iterations to ensure the $\phi$ of the resulting output coordinate matches the $\phi$ used as input:
fs_norm_constant_phi = SBGeom.Flux_Surfaces_Normal_Extended_Constant_Phi(fs)
rzphi = fs_norm.Return_Cylindrical_Position( s = s, LCFS_distance_label= d, theta = 0.2 * np.pi, phi = 0.5 * np.pi)
rzphi_cphi = fs_norm_constant_phi.Return_Cylindrical_Position(s = s, LCFS_distance_label= d, theta = 0.2 * np.pi, phi = 0.5 * np.pi)
print("Normal version phi output - phi input:", rzphi[2] - 0.5 * np.pi)
print("Constant phi version phi output - phi input:", rzphi_cphi[2] - 0.5 * np.pi)
Normal version phi output - phi input: 0.011577579928954584 Constant phi version phi output - phi input: 4.762856775641922e-13
surface_cphi = fs_norm_constant_phi.Mesh_Surface(s, d, N_lines_theta = ntheta, N_lines_phi = nphi, phi_min = phi_min, phi_max = phi_max)
#=================================================================================================
# Plotting:
#==================================================================================================
fig = create_figure_planes()
Mesh.Plot(fig, surface_cphi, color = 'green', name = "Constant phi ")
Mesh.Plot(fig, surface, color = 'red', name = "Normal vector ")
Mesh.Plot(fig, lcfs, color = 'rgb(200, 130, 200)', name= "LCFS", opacity = 0.8)
fig.show()
Clearly, the two approaches define the same surface, but one has the nice property that $\phi_{in} = \phi_{out}$!
Converting into Fourier surfaces¶
In order to have a representation more like the inner flux surfaces, it is also desirable to be able to convert such surfaces to a VMEC Fourier representation.
Internally, the Discrete Fourier Transform is used (via the fft) to obtain the coefficients, which are then converted to the VMEC format. This assumes coordinates equispaced in $\theta ,\phi$
As the normal vector sampling curve is not actually on the $\phi$ input, the resulting surface will not match the desired surface exactly.
This is shown in the below plots: although both Fourier fitted extension surfaces are now exactly $\phi$ slices, the normal fitted does not match the exact surface.
nu_sample = 50
nv_sample = 51
fs_fourier = VMEC.Convert_to_Fourier_Extended(fs_norm, [d], nu_sample, nv_sample)
fs_fourier_cp = VMEC.Convert_to_Fourier_Extended(fs_norm_constant_phi, [d], nu_sample, nv_sample)
surface_cphi_fit = fs_fourier_cp.Mesh_Surface(s, d, N_lines_theta = ntheta, N_lines_phi = nphi, phi_min = phi_min, phi_max = phi_max)
surface_normal_fit = fs_fourier.Mesh_Surface( s, d, N_lines_theta = ntheta, N_lines_phi = nphi, phi_min = phi_min, phi_max = phi_max)
#==================================================================================================================================================================================================================================
# Plotting:
#==================================================================================================================================================================================================================================
fig = create_figure_planes()
Mesh.Plot(fig, surface_cphi_fit, color = 'green', name = "Constant phi ")
Mesh.Plot(fig, surface_normal_fit, color = 'red', name = "Normal vector ")
Mesh.Plot(fig, lcfs, color = 'rgb(200, 130, 200)', name= "LCFS", opacity = 0.8)
fig.show()
def create_figure_equal():
fig = go.Figure()
fig.update_layout(autosize=False, width=600, height=500)
fig.update_xaxes(scaleanchor = "y", scaleratio = 1)
return fig
fig = create_figure_equal()
VMEC.Plot_Poloidal_Slice(fig, fs = fs_norm_constant_phi, s = s, LCFS_distance_label=d, phi=0.0, line_kwargs=dict(color = 'black', width = 4, dash = 'dot'), showlegend=True, name = "Desired Surface")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp, s = s, LCFS_distance_label=d, phi=0.0, line_kwargs=dict(color = 'green', width = 2), showlegend=True, name = "Constant phi Fourier fit")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier, s = s, LCFS_distance_label=d, phi=0.0, line_kwargs=dict(color = 'red', width = 2), showlegend=True, name = "Normal Vector Fit")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_norm_constant_phi, s = s, LCFS_distance_label=d, phi= 2 * np.pi / SYMM / 2, line_kwargs=dict(color = 'black', width = 4, dash = 'dot'))
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp, s = s, LCFS_distance_label=d, phi= 2 * np.pi / SYMM / 2, line_kwargs=dict(color = 'green', width = 2))
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier, s = s, LCFS_distance_label=d, phi= 2 * np.pi / SYMM / 2, line_kwargs=dict(color = 'red', width = 2))
fig.show()
Meshing¶
The resulting Fourier surfaces can be meshed to use for simulations:
fig = go.Figure()
nplasma = 10
nblanket = 10
ntheta = 30
nphi = 30
s_mesh = np.concatenate([np.linspace(0.4,1,nplasma)**2, np.linspace(1,1,nblanket)])
d_mesh = np.concatenate([np.linspace(0.0,0,nplasma) , np.linspace(0,d,nblanket)])
outer_Edge = fs_fourier_cp.Mesh_Surfaces_Closed(s_1 = s_mesh[0], s_2 = s_mesh[-1], LCFS_distance_label_1=d_mesh[0], LCFS_distance_label_2=d_mesh[-1],
N_lines_theta =ntheta, N_lines_phi = nphi, phi_min = phi_min, phi_max = phi_max)
#=================================================================================
# Plotting:
#=================================================================================
fig = create_figure_planes()
Mesh.Plot(fig, outer_Edge, wireframe=True, surface=True, color="red", name="Mesh")
fig.show()
However, this mesh is rather imbalanced: some triangles have sizes much greater htan others, especially near the edge. This is not necessarily a big problem in ray-tracing transport calculations, as you can make the mesh much finer without much cost, however, in deterministic, meshed simulations this is a big issue: convergence is limited by the size, so the mesh needs to be refined. Although, in principle, this could also be done locally, it is much easier to ensure each element is roughly the same size instead.
Equal arc length fitting¶
Since each element is necessarily the same size in the $\phi$ direction (since we define the spacing to be uniform in that direction), and similarly, in the distance away from the LCFS it is controllable, the only direction not uniformly spaced is the $u$ direction.
Luckily, the parametrisation of this coordinate is arbitrary: it doesn't represent a physical coordinate directly like $\phi$. So, we can exploit this parametrisation in the fitting of the Fourier surfaces to obtain a Fourier parametrisation which will have a more even mesh size.
Technically, it works by at each $\phi$ slice (at each sample $\phi$) calculating what $\theta$ values would result in an equally spaced arclength, and then Fourier transforming the sampled points as if they were equally spaced $\theta$ values.
nu_sample = 50
nv_sample = 51
fs_fourier_cp_equal = VMEC.Convert_to_Fourier_Extended(fs_norm_constant_phi, [d], nu_sample, nv_sample, equal_arclength=True)
outer_Edge = fs_fourier_cp_equal.Mesh_Surfaces_Closed(s_1 = s_mesh[0], s_2 = s_mesh[-1], LCFS_distance_label_1=d_mesh[0], LCFS_distance_label_2=d_mesh[-1], N_lines_theta =ntheta, N_lines_phi = nphi, phi_min = phi_min, phi_max = phi_max)
#==================================================================================================================================================================================================================================
# Plotting:
#==================================================================================================================================================================================================================================
fig = create_figure_planes()
Mesh.Plot(fig, outer_Edge, wireframe=True, surface=True, color="blue", name="Equal arclength")
fig.show()
fig = create_figure_equal()
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp, s = s, LCFS_distance_label=d, phi= 0.0, line_kwargs=dict(color = 'red', width = 2 ), nu = ntheta, mode='lines+markers', showlegend=True, name = "Non-equal arclength")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp_equal, s = s, LCFS_distance_label=d, phi= 0.0, line_kwargs=dict(color = 'blue', width = 2 ), nu = ntheta, mode='lines+markers', showlegend=True, name = "Equal arclength")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp, s = s, LCFS_distance_label=d, phi= 2 * np.pi / SYMM / 2, line_kwargs=dict(color = 'red', width = 2 ), mode="lines+markers", nu = ntheta)
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp_equal, s = s, LCFS_distance_label=d, phi= 2 * np.pi / SYMM / 2, line_kwargs=dict(color = 'blue', width = 2 ), mode="lines+markers", nu = ntheta)
fig.show()
As a bonus, the equal arc-length method of fitting has a considerably lower bandwidth in the Fourier space. Truncating the expansion beyond a Fourier harmonic is much less bad for the surface quality in the equal arc length case:
(as a note: if you are intending to truncate in any case, it is better to just sample less than to sample more and truncate)
def Plot_Both_Equal_Non_Equal(fig, ntor, mpol):
nu_sample = 49
nv_sample = 51
fs_fourier_cp_equal_trunc = VMEC.Convert_to_Fourier_Extended(fs_norm_constant_phi, [d], nu_sample, nv_sample, equal_arclength=True, ntor=ntor, mpol=mpol)
fs_fourier_cp_trunc = VMEC.Convert_to_Fourier_Extended(fs_norm_constant_phi, [d], nu_sample, nv_sample, equal_arclength=False, ntor=ntor, mpol=mpol)
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp_equal_trunc, s = s, LCFS_distance_label=d, phi= 2 * np.pi / 4/ SYMM, line_kwargs=dict(color = 'red', width = 2 ), nu = ntheta, mode='lines+markers', showlegend=True, name = "Equal arclength")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_fourier_cp_trunc, s = s, LCFS_distance_label=d, phi= 2 * np.pi / 4/ SYMM, line_kwargs=dict(color = 'blue', width = 2 ), nu = ntheta, mode='lines+markers', showlegend=True, name = "Non-equal arclength")
VMEC.Plot_Poloidal_Slice(fig, fs = fs_norm_constant_phi, s = s, LCFS_distance_label=d, phi= 2 * np.pi / 4/ SYMM, line_kwargs=dict(color = 'black', width = 2, dash='dot' ), nu = 300, mode='lines', showlegend=True, name = "Desired Curve")
fig = create_figure_equal()
Plot_Both_Equal_Non_Equal(fig, None, None )
fig.show()
fig = create_figure_equal()
Plot_Both_Equal_Non_Equal(fig, 6, 6)
fig.show()
Non-Uniform distance¶
There are also tools supporting the use of non-uniform distances.
To use this, one would have to define a thickness matrix, having the distances at a set of toroidal and poloidal coordinates.
tor_max = 2 * np.pi / SYMM / 2
toroidal_angles = np.array([0.0, tor_max / 4, tor_max / 2, tor_max / 4 * 3, tor_max])
poloidal_angles = np.array([0.0, np.pi / 2, np.pi, 3 * np.pi / 2 , 2 * np.pi])
thickness_matrix =np.array([
[1.4, 1.0, 0.6, 1.0, 1.4], # toroidal 0
[1.2, 0.8, 0.6, 0.8, 1.2], # toroidal max / 4
[1.0, 0.7, 0.6, 0.7, 1.0], # toroidal max / 2
[1.2, 0.8, 0.6, 0.8, 1.2], # toroidal max * 3 / 4
[0.6, 0.6, 0.6, 0.6, 0.6] # toroidal max
]).T # transpose because poloidal should be first coordinate (row of this matrix)
fs_ext = VMEC.Convert_to_Fourier_Extended(fs_norm_constant_phi, [1.0], 100 , 101, LCFS_distance_labels_functions=[VMEC.Fit_From_Thickness_Matrix(poloidal_angles, toroidal_angles, thickness_matrix, method='cubic')],equal_arclength=True, ntor = None, mpol = None)
#==================================================================================================================================================================================================================================
# Plotting:
#==================================================================================================================================================================================================================================
nu = 50
nv = 51
ext_surf = fs_ext.Mesh_Surface(1.0,1.0,nu,nv, 0.0 , tor_max )
surf = fs_ext.Mesh_Surface(1.0,0.0,50,60, -0.05, tor_max + 0.05)
fig = create_figure_planes()
Mesh.Plot(fig, ext_surf, wireframe = True, name = "Non-uniform distance", showlegend=True)
Mesh.Plot(fig, surf , name = "LCFS", showlegend=True, color="purple", opacity =0.6)
fig.show()
fig = create_figure_equal()
phi = tor_max
def plot_phi(phi, showlegend=False):
VMEC.Plot_Poloidal_Slice(fig, fs_norm_constant_phi, 1.0, 0.6, phi,line_kwargs = dict(color = 'red'), name='Uniform d=0.6', showlegend=showlegend)
VMEC.Plot_Poloidal_Slice(fig, fs_ext, 1.0, 1.0, phi,line_kwargs = dict(color = 'lightblue'), name ="Non-uniform distance", showlegend=showlegend)
VMEC.Plot_Poloidal_Slice(fig, fs_ext, 1.0, 0.0, phi,line_kwargs = dict(color = 'purple'), name ="LCFS", showlegend=showlegend)
plot_phi(tor_max, showlegend=True)
plot_phi(0.1)
plot_phi(tor_max * 0.563)
fig.show()