In this example we demonstrate how to compute a current potential on one surface that minimizes the error of the normal magnetic field on another. This is relevant to stellarator optimization when the normal field is minimized on a desired plasma surface and the current potential lies on a so-called winding surface here coils will be initialized for further optimization.

Note: Running this example requires the h5py python package, which is installable using pip or other standard methods.

import struct
import sys
import os
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
plt.rcParams['figure.figsize']=(6,6)
plt.rcParams['font.weight']='bold'
plt.rcParams['axes.labelweight']='bold'
plt.rcParams['lines.linewidth']=2
plt.rcParams['lines.markeredgewidth']=2
%matplotlib inline
%config InlineBackend.figure_format = "retina"

Load ThinCurr library

To load the ThinCurr python module we need to tell python where to the module is located. This can be done either through the PYTHONPATH environment variable or within a script using sys.path.append() as below, where we look for the environement variable OFT_ROOTPATH to provide the path to where the OpenFUSIONToolkit is installed (/Applications/OFT for binaries on macOS).

thincurr_python_path = os.getenv('OFT_ROOTPATH')
if thincurr_python_path is not None:
    sys.path.append(os.path.join(thincurr_python_path,'python'))
from OpenFUSIONToolkit.ThinCurr import ThinCurr
from OpenFUSIONToolkit.ThinCurr.meshing import write_ThinCurr_mesh, build_regcoil_grid, build_periodic_mesh, write_periodic_mesh
from OpenFUSIONToolkit.util import build_XDMF

Generate mesh and ThinCurr input for coil and plasma surfaces for NCSX

First we create ThinCurr models for the winding and plasma surfaces using a REGCOIL definition of the NCSX stellarator using build_regcoil_grid(). This subroutine builds a uniform grid over one field period, which can then be used by build_periodic_mesh() to build a mesh from the resulting grid. The result is a ThinCurr model, including periodicity mapping information (r_map).

# Create coil mesh
coil_grid, nfp = build_regcoil_grid('regcoil_out.li383.nc', 'coil', 64, 64, full_torus=False)
lc_coil, r_coil, tnodeset_coil, pnodesets_coil, r_map = build_periodic_mesh(coil_grid,nfp)
# Create plasma mesh
plasma_grid, _ = build_regcoil_grid('regcoil_out.li383.nc', 'plasma', 64, 64*nfp, full_torus=True)
lc_plasma, r_plasma, tnodeset_plasma, [pnodeset_plasma], _ = build_periodic_mesh(plasma_grid,1)

Plot resulting winding surface mesh

fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1, projection='3d')
ax.plot(r_coil[tnodeset_coil,0], r_coil[tnodeset_coil,1], r_coil[tnodeset_coil,2], c='red')
for pnodeset_coil in pnodesets_coil:
    ax.plot(r_coil[pnodeset_coil,0], r_coil[pnodeset_coil,1], r_coil[pnodeset_coil,2], c='blue')
 
ax = fig.add_subplot(1,2,2, projection='3d')
_ = ax.plot_trisurf(r_coil[:,0], r_coil[:,1], r_coil[:,2], triangles=lc_coil, cmap='viridis')

Plot resulting plasma surface mesh

fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1, projection='3d')
ax.plot(r_plasma[tnodeset_plasma,0], r_plasma[tnodeset_plasma,1], r_plasma[tnodeset_plasma,2], c='red')
ax.plot(r_plasma[pnodeset_plasma,0], r_plasma[pnodeset_plasma,1], r_plasma[pnodeset_plasma,2], c='blue')
 
ax = fig.add_subplot(1,2,2, projection='3d')
_ = ax.plot_trisurf(r_plasma[:,0], r_plasma[:,1], r_plasma[:,2], triangles=lc_plasma, cmap='viridis')

Save models to HDF5 format

write_ThinCurr_mesh('thincurr_plasma.h5', r_plasma, lc_plasma+1, np.ones((lc_plasma.shape[0],)),
                    holes=[pnodeset_plasma+1,tnodeset_plasma+1])
write_periodic_mesh('thincurr_coil.h5', r_coil, lc_coil+1, np.ones((lc_coil.shape[0],)),
                    tnodeset_coil, pnodesets_coil, pmap=r_map, nfp=nfp, include_closures=False)

Saving mesh: thincurr_plasma.h5

Saving mesh: thincurr_coil.h5

Compute optimal winding surface current

Setup ThinCurr model for winding surface

We now create a ThinCurr instance for the winding surafce. Once created, we setup the model from the existing HDF5 mesh definition generated above using setup_model(). We also initialize I/O for this model using setup_io() to enable output of plotting files for 3D visualization in VisIt or Paraview.

In this case we specify a coil directory to use for saving I/O files to keep things separate for the other cases to be run in this notebook and in ThinCurr Python Example: Compute frequency-response in a torus.

Note: A warning will be generated that no XML node was found and a default resistivity value is being used. This is OK as we will not use the resitivity matrix with this model.

tw_coil = ThinCurr(nthreads=4)
tw_coil.setup_model(mesh_file='thincurr_coil.h5')
tw_coil.setup_io(basepath='coil')

#----------------------------------------------
Open FUSION Toolkit Initialized
Development branch:   main
Revision id:          8440e61
Parallelization Info:
  Not compiled with MPI
  # of OpenMP threads =    4
Fortran input file    = oftpyin                                                                         
XML input file        = none                                                                            
Integer Precisions    =    4   8
Float Precisions      =    4   8  16
Complex Precisions    =    4   8
LA backend            = native
#----------------------------------------------


Creating thin-wall model
 Orientation depth =       24192
  Loading V(t) driver coils
  Loading I(t) driver coils

  # of points    =        12096
  # of edges     =        36288
  # of cells     =        24192
  # of holes     =            4
  # of Vcoils    =            0
  # of closures  =            0
  # of Icoils    =            0

  Building holes


WARNING: Unable to find "thincurr" XML node
WARNING: No "thincurr" XML node, using "eta=mu0" for all regions

Setup ThinCurr model for plasma surface

We now create a ThinCurr instance for the winding surafce. Once created, we setup the model from the existing HDF5 mesh definition generated above using setup_model(). We also initialize I/O for this model using setup_io() to enable output of plotting files for 3D visualization in VisIt or Paraview.

In this case we specify a plasma directory to use for saving I/O files to keep things separate for the other cases to be run in this notebook and in ThinCurr Python Example: Compute frequency-response in a torus.

Note: A warning will be generated that no XML node was found and a default resistivity value is being used. This is OK as we will not use the resitivity matrix with this model.

tw_plasma = ThinCurr()
tw_plasma.setup_model(mesh_file='thincurr_plasma.h5')
tw_plasma.setup_io(basepath='plasma')

Creating thin-wall model
 Orientation depth =       24576
  Loading V(t) driver coils
  Loading I(t) driver coils

  # of points    =        12288
  # of edges     =        36864
  # of cells     =        24576
  # of holes     =            2
  # of Vcoils    =            0
  # of closures  =            0
  # of Icoils    =            0

  Building holes


WARNING: Unable to find "thincurr" XML node
WARNING: No "thincurr" XML node, using "eta=mu0" for all regions

Compute mutual inductance between coil and plasma surfaces

coupling = tw_coil.cross_coupling(tw_plasma,cache_file='Lmat_coupling.save')
nelems_new = coupling.shape[0]-nfp+1
coupling_new = np.zeros((nelems_new,coupling.shape[1]))
coupling_new[:-1,:] = coupling[:-nfp,:]
coupling_new[-1,:] = coupling[-nfp:,:].sum(axis=0)
hole_start = coupling_new.shape[1]-2
del coupling

 Building element<->element mutual inductance matrix
     Time =  1m  9s

Compute coil current regularization matrix

reg_mat = tw_coil.get_regmat()
ncold = int(reg_mat.shape[1]/nfp)
reg_new = np.zeros((nelems_new,ncold))
reg_new[:-1,:] = reg_mat[:-nfp,:ncold]
reg_new[-1,:] = reg_mat[-nfp:,:ncold].sum(axis=0)
print(reg_new[0,:].max())
del reg_mat

1.2828298118544255

Compute optimal current distribution with specified regularization

lam = 1.E-4
wt1 = 1.E-13
wt2 = 2.E-3
#
reg_coupling = np.hstack((coupling_new,lam*reg_new))
lhs = np.zeros((reg_coupling.shape[1],))
#
lhs[hole_start] = 1.E0
reg_coupling[:,hole_start] *= wt1/lhs[hole_start]
lhs[hole_start] *= wt1/lhs[hole_start]
#
lhs[hole_start+1] = 1.E0
reg_coupling[:,hole_start+1] *= wt2/lhs[hole_start+1]
lhs[hole_start+1] *= wt2/lhs[hole_start+1]
#
reg_coupling[:,-2]=0.0; reg_coupling[-2,-2]=1.0
#
ls_output = np.linalg.lstsq(np.transpose(reg_coupling), lhs, rcond=None)
rhs = ls_output[0]

Save error and current distribution for plotting

# Save error
tmp_err = np.dot(np.transpose(coupling_new),rhs)
tw_plasma.save_scalar(tmp_err[:-2],'E')
build_XDMF(path='plasma')
# Expand current potential to holes and save
output = np.zeros((nelems_new+nfp-1,))
output[:nelems_new] = rhs
output[-nfp:] = rhs[-1]
tw_coil.save_current(output,'J')
build_XDMF(path='coil')

Removing old Xdmf files
Creating output files
Removing old Xdmf files
Creating output files

Plot current potential on one field period

data_unwrap=np.r_[0.0,rhs[:-2]].reshape((63,64))
plt.contour(data_unwrap.transpose(),20)
_ = plt.colorbar()

Plot normali field error on one field period

data_unwrap=tmp_err[:64*64].reshape((64,64))
plt.contourf(data_unwrap.transpose(),20)
_ = plt.colorbar()