In this example we show how to compute equilibria in ITER with H-mode like profiles for the "inverse" case where we have a desired shape, plasma current and pressure, but the required coil currents are unkown.

This example utilizes the mesh built in TokaMaker Meshing Example: Building a mesh for ITER.

Note: Running this example requires the h5py python packages. Both of which are installable using pip or other standard methods. This example also requires use of the Sauter formula for bootstrap current (See: O. Sauter, et al., Phys. Plasmas 6, 2834 (1999); DOI:10.1063/1.873240). It is available as part of the open source platform OMFIT, written as part of the Kolemen Group Automatic Kinetic EFIT Project "auto_kEFIT" (See: https://omfit.io/_modules/omfit_classes/utils_fusion.html).

This formula can be installed using pip: $ pip install --upgrade omfit_classes

The utility function Hmode_profiles() is also imported from omfit_classes to facilitate specification of a pressure pedestal.

import os
import sys
import numpy as np
import matplotlib.pyplot as plt
from omfit_classes.utils_fusion import sauter_bootstrap, Hmode_profiles
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 TokaMaker library

To load the TokaMaker python module we need to tell python where to the module is located. This can be done either through the PYTHONPATH environment variable or using 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 on macOS).

tokamaker_python_path = os.getenv('OFT_ROOTPATH')
if tokamaker_python_path is not None:
    sys.path.append(os.path.join(tokamaker_python_path,'python'))
from OpenFUSIONToolkit import OFT_env
from OpenFUSIONToolkit.TokaMaker import TokaMaker, solve_with_bootstrap
from OpenFUSIONToolkit.TokaMaker.meshing import load_gs_mesh

Compute equilibria

Initialize TokaMaker object

We now create a OFT_env instance for execution using two threads and a TokaMaker instance to use for equilibrium calculations. Note at present only a single TokaMaker instance can be used per python kernel, so this command should only be called once in a given Jupyter notebook or python script. In the future this restriction may be relaxed.

myOFT = OFT_env(nthreads=2)

mygs = TokaMaker(myOFT)

#----------------------------------------------
Open FUSION Toolkit Initialized
Development branch:   v1_beta6
Revision id:          681e857
Parallelization Info:
  # of MPI tasks      =    1
  # of NUMA nodes     =    1
  # of OpenMP threads =    2
Fortran input file    = /var/folders/52/n5qxh27n4w19qxzqygz2btbw0000gn/T/oft_64680/oftpyin
XML input file        = none
Integer Precisions    =    4   8
Float Precisions      =    4   8  16
Complex Precisions    =    4   8
LA backend            = native
#----------------------------------------------

Load mesh into TokaMaker

Now we load the mesh generated in TokaMaker Meshing Example: Building a mesh for ITER using load_gs_mesh() and setup_mesh. Then we use setup_regions(), passing the conductor and coil dictionaries for the mesh, to define the different region types. Finally, we call setup() to setup the required solver objects. During this call we can specify the desired element order (min=2, max=4) and the toroidal field through F0 = B0*R0, where B0 is the toroidal field at a reference location R0.

mesh_pts,mesh_lc,mesh_reg,coil_dict,cond_dict = load_gs_mesh('ITER_mesh.h5')
mygs.setup_mesh(mesh_pts, mesh_lc, mesh_reg)
mygs.setup_regions(cond_dict=cond_dict,coil_dict=coil_dict)
mygs.settings.maxits=100
mygs.setup(order = 2, F0 = 5.3*6.2)

**** Loading OFT surface mesh

**** Generating surface grid level  1
  Generating boundary domain linkage
  Mesh statistics:
    Area         =  2.859E+02
    # of points  =    4757
    # of edges   =   14156
    # of cells   =    9400
    # of boundary points =     112
    # of boundary edges  =     112
    # of boundary cells  =     112
  Resolution statistics:
    hmin =  9.924E-03
    hrms =  2.826E-01
    hmax =  8.466E-01
  Surface grounded at vertex     870


**** Creating Lagrange FE space
  Order  =    2
  Minlev =   -1

 Computing flux BC matrix 
 Inverting real matrix
   Time =    1.1390000000000000E-003

Define a vertical stability coil

Like many elongated equilibria, the equilibrium we seek to compute below is vertically unstable. In this case we use the actual ITER Vertical Stability Coil (VSC) in order to help with convergence using the set_coil_vsc() method.

Note: While ITER has a "real" VSC, this is not required and this functionality can instead be used to define a "virtual" VSC by pairing coils in a way that are not necessarily paired experimentally.

mygs.set_coil_vsc({'VS': 1.0})

Define hard limits on coil currents

Hard limits on coil currents can be set using set_coil_bounds(). In this case we just the simple and approximate bi-directional limit of 50 MA in each coil.

Bounds are specified using a dictionary of 2 element lists, containing the minimum and maximum bound, where the dictionary key corresponds to the coil names, which are available in mygs.coil_sets

coil_bounds = {key: [-50.E6, 50.E6] for key in mygs.coil_sets}

mygs.set_coil_bounds(coil_bounds)

Compute Inverse Equilibrium

Define global quantities and targets

For the inverse case we define a target for the plasma current and the peak plasma pressure, which occurs on the magnetic axis.

Note: These constraints can be considered "hard" constraints, where they will be matched to good tolerance as long as the calculation converges.

Ip_target=15.2E6
P0_target=6.2E5
mygs.set_targets(Ip=Ip_target, pax=P0_target)

Define shape targets

In order to constrain the shape of the plasma we can utilize two types of constraints:

isoflux points, which are points we want to lie on the same flux surface (eg. the LCFS)
saddle points, where we want the poloidal magnetic field to vanish (eg. X-points). While one can also use this constraint to enforce a magnetic axis location, instead set_targets() should be used with arguments R0 and V0.

Note: These constraints can be considered "soft" constraints, where the calculation attempts to minimize error in satisfying these constraints subject to other constraints and regularization.

Here we define a handful of isoflux points that we want to lie on the LCFS of the target equilibrium. Additionally, we define a single X-point and set it as a saddle constraint as well as adding it to the list of isoflux points.

isoflux_pts = np.array([
    [ 8.20,  0.41],
    [ 8.06,  1.46],
    [ 7.51,  2.62],
    [ 6.14,  3.78],
    [ 4.51,  3.02],
    [ 4.26,  1.33],
    [ 4.28,  0.08],
    [ 4.49, -1.34],
    [ 7.28, -1.89],
    [ 8.00, -0.68]
])
x_point = np.array([[5.125, -3.4],])
mygs.set_isoflux(np.vstack((isoflux_pts,x_point)))
mygs.set_saddles(x_point)

Define coil regularization matrix

In general, for a given coil set a given plasma shape cannot be exactly reproduced, which generally yields large amplitude coil currents if no constraint on the coil currents is applied. As a result, it is useful to include regularization terms for the coils to balance minimization of the shape error with the amplitude of current in the coils. In TokaMaker these regularization terms have the general form, where each term corresponds to a set of coil coefficients, target value, and weight. The coil_reg_term() method is provided to aid in defining these terms.

In this case, one regularization term is added for each coil with a single unit coefficient for that coil and target of zero with modest weights. This regularization acts to penalize the amplitude of current in each coil, acting to balance coil current with error in the shape targets. Note that the weight on the VSC virtual coil (#VSC) defined above is set high to prevent interaction with the real VS coil set (see below for further information).

# Set regularization weights
regularization_terms = []
for name, coil in mygs.coil_sets.items():
    # Set zero target current and different small weights to help conditioning of fit
    if name.startswith('CS'):
        if name.startswith('CS1'):
            regularization_terms.append(mygs.coil_reg_term({name: 1.0},target=0.0,weight=2.E-2))
        else:
            regularization_terms.append(mygs.coil_reg_term({name: 1.0},target=0.0,weight=1.E-2))
    elif name.startswith('PF'):
        regularization_terms.append(mygs.coil_reg_term({name: 1.0},target=0.0,weight=1.E-2))
    elif name.startswith('VS'):
        regularization_terms.append(mygs.coil_reg_term({name: 1.0},target=0.0,weight=1.E-2))
# Disable VSC virtual coil
regularization_terms.append(mygs.coil_reg_term({'#VSC': 1.0},target=0.0,weight=1.E2))
 
# Pass regularization terms to TokaMaker
mygs.set_coil_reg(reg_terms=regularization_terms)

Define flux functions

Although TokaMaker has a "default" profile for the F*F' and P' terms this should almost never be used and one should instead choose an appropriate flux function for their application. In this case we use an H-mode-like kinetic profiles parametrized via the OMFIT function Hmode_profiles() (see note at top). Within TokaMaker these profiles are represented as piecewise linear functions, which can be set up using the dictionary approach shown below.

Using the OMFIT Hmode_profiles() profile parametrization function, we can specify ITER-like H-mode profiles. See Figs. 18 and 19 for these $n_e$ and $T_e$ profiles: https://iopscience.iop.org/article/10.1088/0029-5515/48/7/075005/pdf

For pressure, we are assuming a quasineutral and isothermal plasma, (i.e. ne = ni, Te = Ti), however this assumption is not usually valid in the core and separate ion profiles will need to be specified.

Note: More tools to aide in setting these profiles are coming soon

n_sample = 100
psi_sample = np.linspace(0.0,1.0,n_sample)
 
 
xphalf = 0.965
widthp_Te = 0.1
widthp_ne = 0.35
 
ne = Hmode_profiles(edge=0.35, ped=0.6, core=1.1, rgrid=n_sample, expin=1.6, expout=1.6, widthp=widthp_ne, xphalf=xphalf) * 1e20
Te = Hmode_profiles(edge=1500., ped=5000., core=21000., rgrid=n_sample, expin=1.3, expout=1.7, widthp=widthp_Te, xphalf=xphalf)
ni = ne.copy() # Assuming quasineutrality
Ti = Te.copy() # Assuming isothermal
pressure = (1.602e-19 * ne * Te) + (1.602e-19 * ni * Ti) # 1.602e-19 * [m^-3] * [eV] = [Pa]

fig, ax = plt.subplots(3,1,sharex=True)
# Plot ne
ax[0].plot(psi_sample,ne)
ax[0].set_ylabel(r'ne [m$^{-3}$]')
ax[0].set_ylim(0.,1.2e+20)
# Plot Te
ax[1].plot(psi_sample,Te)
ax[1].set_ylabel("Te [eV]")
ax[1].set_ylim(0.,23000)
# Plot total pressure
ax[2].plot(psi_sample,pressure)
ax[2].set_ylabel("Pressure [Pa]")
ax[-1].set_xlabel(r"$\hat{\psi}$")
_ = ax[-1].set_xlabel(r"$\hat{\psi}$")

Using the OMFIT Hmode_profiles() profile parametrization function, we can specify an inductive, L-mode-like $j_{tor}$ profile. See Fig. 7 for a reference $j_{total}$ profile including the bootstrap contribution: https://iopscience.iop.org/article/10.1088/0029-5515/48/7/075005/pdf

This $j_{inductive}$ profile was generated by hand such that the sum of $j_{inductive}$ and the boostrap current profile ( $j_{bootstrap}$) calculated by the Sauter formula approximately equal the high fidelity $j_{total}$ profile in Fig. 7. Achieving this may require referencing a published $j_{total}$ profile which was calculated via transport simulations in similar conditions to your goal TokaMaker equilibrium, some iteration by hand, or a combination of the two.

jtor_noBS = Hmode_profiles(edge=0.0, ped=0.001, core=0.7, rgrid=n_sample, expin=0.9, expout=3.5, widthp=0.001, xphalf=None) * 1e+6 # [A/m^2]

fig, ax = plt.subplots(1,1,sharex=True)
# Plot j_inductive
ax.plot(psi_sample,jtor_noBS,label='Inductive $j_{tor}$')
ax.set_ylabel("j [A/m$^2$]")
ax.set_xlabel(r"$\hat{\psi}$")
ax.set_title("H-mode $j_{inductive}$")
_ = ax.legend()

Compute equilibrium

We can now compute a free-boundary equilibrium using these constraints. Note that before running a calculation for the first time we must initialize the flux function $\psi$, which can be done using init_psi(). This subroutine initializes the flux using the specified Ip_target from above, which is evenly distributed over the entire plasma region or only with a boundary defined using a center point (R,Z), minor radius (a), and elongation and triangularity. Coil currents are also initialized at this point using the constraints above and this uniform plasma current initialization.

solve() is then called to compute a self-consitent Grad-Shafranov equilibrium. If the result variable (err_flag) is zero then the solution has converged to the desired tolerance ( $10^{-6}$ by default).

R0 = 6.3
Z0 = 0.5
a = 2.0
kappa = 1.4
delta = 0.0
err_flag = mygs.init_psi(R0, Z0, a, kappa, delta)
err_flag = mygs.solve()

Starting non-linear GS solver
     1  1.5142E+00  9.9086E-02  5.9853E-01  6.4538E+00  5.3436E-01  5.3525E-04
     2  4.6368E+00  6.0157E-02  1.5071E-01  6.4193E+00  5.3221E-01  6.8799E-04
     3  4.9253E+00  5.6147E-02  3.6043E-02  6.4024E+00  5.3113E-01  7.0570E-04
     4  4.9593E+00  5.5547E-02  8.9613E-03  6.3957E+00  5.3065E-01  7.0442E-04
     5  4.9593E+00  5.5481E-02  2.5544E-03  6.3933E+00  5.3036E-01  7.0151E-04
     6  4.9569E+00  5.5487E-02  9.1473E-04  6.3925E+00  5.3019E-01  6.9957E-04
     7  4.9556E+00  5.5495E-02  3.8318E-04  6.3922E+00  5.3011E-01  6.9853E-04
     8  4.9551E+00  5.5499E-02  1.6704E-04  6.3921E+00  5.3007E-01  6.9803E-04
     9  4.9549E+00  5.5501E-02  7.2262E-05  6.3921E+00  5.3004E-01  6.9780E-04
    10  4.9548E+00  5.5502E-02  3.0753E-05  6.3921E+00  5.3004E-01  6.9770E-04
    11  4.9548E+00  5.5502E-02  1.2905E-05  6.3921E+00  5.3003E-01  6.9766E-04
    12  4.9548E+00  5.5502E-02  5.3627E-06  6.3921E+00  5.3003E-01  6.9764E-04
    13  4.9548E+00  5.5502E-02  2.2167E-06  6.3921E+00  5.3003E-01  6.9763E-04
    14  4.9548E+00  5.5502E-02  9.1477E-07  6.3921E+00  5.3003E-01  6.9763E-04
 Timing:  0.10165199998300523     
   Source:     3.4617999859619886E-002
   Solve:      3.6162999866064638E-002
   Boundary:   3.4109998960047960E-003
   Other:      2.7460000361315906E-002

Now input kinetic profiles to obtain H-mode solution

solve_with_bootstrap() is now called to compute a self-consistent Grad-Shafranov equilibrium that includes non-inductive current (the bootstrap current). Inclusion of bootstrap current effects is essential for H-mode equilibrium stability. If the result variable (err_flag) is zero then the solution has converged to the desired tolerance ( $10^{-6}$ by default).

Zeff=np.linspace(1.5,1.5,n_sample)
ni = ne.copy()
Ti = Te.copy()
 
mygs.settings.pm=False
mygs.update_settings()
mygs.set_targets(pax=pressure[0])
 
err_flag, j_BS = solve_with_bootstrap(mygs,ne,Te,ni,Ti,jtor_noBS,Zeff)

>>> Initializing equilibrium with pedestal removed:
  Solve flag:  -1
>>> Iterating on H-mode equilibrium solution:
> Iteration 0:
  Solve flag:  0
> Iteration 1:
  Solve flag:  0
> Iteration 2:
  Solve flag:  0

Plot equilibrium

Flux surfaces of the computed equilibrium can be plotted using the plot_psi() method. The additional plotting methods plot_machine() and plot_constraints() are also used to show context and other information. Each method has a large number of optional arguments for formatting and other options.

fig, ax = plt.subplots(1,1)
mygs.plot_machine(fig,ax,coil_colormap='seismic',coil_symmap=True,coil_scale=1.E-6,coil_clabel=r'$I_C$ [MA]')
mygs.plot_psi(fig,ax,xpoint_color=None,vacuum_nlevels=4)
mygs.plot_constraints(fig,ax,isoflux_color='tab:red',isoflux_marker='o')

Plot final self-consistent P' and FF' including $j_{bootstrap}$ contribution

psi,f,fp,p,pp = mygs.get_profiles(npsi=n_sample)
 
fig, ax = plt.subplots(2,1,sharex=True)
# Plot F*F'
ax[0].plot(psi,f*fp)
ax[0].set_ylabel("FF'")
# Plot P'
ax[1].plot(psi,pp)
ax[1].set_ylabel("P'")
_ = ax[-1].set_xlabel(r"$\hat{\psi}$")
 
ax[0].set_title("H-mode P' and FF'")

Text(0.5, 1.0, "H-mode P' and FF'")

Plot final $j_{total}$ = $j_{inductive}$ + $j_{bootstrap}$

fig, ax = plt.subplots(1,1,sharex=True)
 
 
_,_,ravgs,_,_,_ = mygs.get_q(npsi=n_sample) # get flux averaged R and 1/R from equilibrium solution
R_avg = ravgs[0] # R
one_over_R_avg = ravgs[1] # 1/R
mu0 = np.pi*4.E-7
tkmkr_jtor = R_avg * pp + one_over_R_avg*(f*fp) / mu0 # Jtor = R*P' + FF' / (mu0*R)
    
# Plot j_total
ax.plot(psi,tkmkr_jtor/1e+6,label='Output $j_{tor}$') # in MA/m^2
# Plot j_bootstrap
ax.plot(psi,j_BS/1e+6,label='$j_{bootstrap}$') # in MA/m^2
ax.set_ylabel("j [MA/m$^2$]")
ax.set_xlabel(r"$\hat{\psi}$")
ax.set_title("H-mode $j_{tor}$")
_ = ax.legend()

Print equilibrium information and coil currents

Basic parameters can be displayed using the print_info() method. For access to these quantities as variables instead the get_stats() can be used.

The final coil currents can also be retrieved using the get_coil_currents() method, which are all within the approximate coil limits imposed above.

mygs.print_info()
 
print()
print("Coil Currents [MA]:")
coil_currents, _ = mygs.get_coil_currents()
for key, current in coil_currents.items():
    print('  {0:10} {1:10.2F}'.format(key+":",current/1.E6))

Equilibrium Statistics:
  Topology                =   Diverted
  Toroidal Current [A]    =    9.0600E+06
  Current Centroid [m]    =    6.452  0.480
  Magnetic Axis [m]       =    6.612  0.520
  Elongation              =    1.870 (U:  1.713, L:  2.027)
  Triangularity           =    0.450 (U:  0.344, L:  0.556)
  Plasma Volume [m^3]     =   803.114
  q_0, q_95               =    2.079  5.393
  Plasma Pressure [Pa]    =   Axis:  7.3989E+05, Peak:  7.3989E+05
  Stored Energy [J]       =    3.3408E+08
  <Beta_pol> [%]          =   177.8486
  <Beta_tor> [%]          =    2.5021
  <Beta_n>   [%]          =    2.8709
  Diamagnetic flux [Wb]   =   -4.6822E-01
  Toroidal flux [Wb]      =    1.1689E+02
  l_i                     =    0.7981

Coil Currents [MA]:
  CS3U:            4.45
  CS2U:            1.72
  CS1U:           -2.37
  CS1L:           -1.71
  CS2L:           -3.80
  CS3L:            4.99
  PF1:             5.46
  PF2:            -0.57
  PF3:            -4.64
  PF4:            -2.59
  PF5:            -7.51
  PF6:            18.55
  VS:             -0.13