PyISOLVER

— A Fast Python OOP Implementation of LRDFIT Algorithm

Version: Beta 0.20

Fanghao Yang, PPPL 10/05/2018

Overview

PyISOVLER is the Python version of ISOLVER ported from the original IDL version of ISOLVER, which is mostly developed by Jon Menard and Bill Davis.

PyISOVLER is not ported from IDL version line by line but totally refactored into OOP style Python code. It contains most frequently used methods and data structures from EFIT, LRDFIT and a few other IDL libraries using by ISOLVER.

Current PyISOVLER is still under development and testing. It only runs with current input data and parameters which are hard coded in the “isolver_input.py” file. It may not run at all or not run correctly for other input data sets.

PyISOVLER is developed with Python 2. However, the code may be moved to Python 3 in near future, Thus, it was written to get best compatibility with Python 3 (does not mean working with Python 3). In this way, only a small effort was need to upgrade to Python 3.

How to Acquire PyISOLVER

To run the PyISOLVER on your local domain, you need use subversion to checkout via:

svn checkout svn+ssh://svnsrv.pppl.gov/svn/nanalysis/pyIsolver/release/beta0.20 ./beta0.20

If you want to change the original code published by Fanghao Yang. You may create a side-line development version in the branch. If necessary, your revision would be merged into the trunk for future release.

How to Run

Now to use it, you need load the module first,

module load nstx/pyisolver

Then, go to the directory of current copy and run the program,

python -s run.py

The maximum number of iterations is limited to 51 for current input parameter. You shall see 5 out of 6 profiles converged within 51 iterations.

How to Change Input Parameters

The “input_generator.py” is a tool for users to generate input JSON files for batch jobs.

For example, if you want to test the convergence for different input sets, you may need change the max iteration number in input JSON file. There are three methods to change input paramters:

1. Edit the function “generate_input”

“generate_input” is a function which has hard-coded input parameters and it will return an instance of Inputs class. You may edit those hard-coded parameters to generate input file. Find where the instance of “IterConvParam” is created in “generate_input()” and “nitmax” is the number to set the maximum iterations.

2. Load existing input file, edit and save the file

You may also load existing JSON input file. For example, you want to load an input parameter file for NSTX.

filename = "input/nstx_test.json"
input_obj = lrdfit.Inputs()
input_obj.load_json(filename)

Then, for example, you want to change the maximum iteration number to 100.

input_obj.itconvpars.nitmax = 100

After you have changed the parameter, you may save the input parameter as JSON file for running ISOLVER code.

input_obj.save_json(filename)

You are free to edit the “input_generator.py” to generate a batch of input files for your job.

3. Edit JSON files

You can also edit the JSON file directly and save it for your job. Since there is no comment in JSON file, you need know what you are changing.

How to Plot Figures

In the example run code “run.py”, in line 40:

ctimes = s.iterate_state_to_convergence(c, m, grs, noplot=False)

If argument “noplot” is false, it will dynamically plot results after each iteration

If argument “noplot” is true, it will disable the plot then no figure would be plotted.

If you want to get faster calculation, then, disable the plot option since plotting takes extra time.

Reference to Original IDL ISOLVER

Current input data sets and profile parameters are derived from the IDL version isolver batch file “input”. For more information, you may find it from the IDL version isolver subversion repository

svn+ssh://svnsrv.pppl.gov/usr/local/svn/isolver/trunk

Example Run

An example of running ISOLVER is included in “run.py”.

The first step is load user input file as JSON format.

First, create an instance of Inputs from lrdfit library. The input data set includes all parameters to initialize ISOLVER and iteratively approach to the solution.

input_obj = lrdfit.Inputs()

Then, use “load_json” method to load the JSON file which stores input data.

input_obj.load_json(filename)

If you have changed some parameters in Inputs data and you may save it as a local JSON file for future work.

input_obj.save_json(filename)

Next, create an instance of GridStructure from lrdfit library for loading and computing grid data with Green’s function.

grs = lrdfit.GridStructure()

Then, use “get_greens_data” method to load the grid data in NetCDF4 format.

grs.get_greens_data(grsfile)

Next, create an instance of FluxCoord from lrdfit library, which has all flux coordinate data for the tokamak device. Some parameters from the Inputs would be used to initialize flux coordinates.

cfc = lrdfit.FluxCoord(nradius=input_obj.fcspars.nradius,
                       ntheta=input_obj.fcspars.ntheta)

The State data comprises all profiles about targeting shot including temperature, pressure, voltage, current, flux and more. Then, an instance of State would be created using Inputs data, moreover, using all existing data sets, the State data would be initialized from Equilibrium data.

The Equilibrium data was previously calculated by EFIT and usually stored as G-EQDSK format. There are two ways to load EFIT G-EQDSK data, either from the local file or from MDS+ tree. All paramters about loading Q-EQDSK data are included as Inputs parameter.

isj = lrdfit.State(inputobj=input_obj)
isj.initialize_from_equilibrium(grs, cfc)

The following step is to calculate the momentum from initial State and flux cooridnate data. The momentum data has similar data structure as class Equilibrium. The momentum data contains boundary data and other data derived from current State profiles and current flux coordinate data.

m = lrdfit.Momentum()
m.compute_state_moments(s, cfc)

The initial State is renamed as State for following iterative computation.

s = isj

Then, the solver starts to solve Grad–Shafranov equation and converge into an acceptable solution. All tolerance of profiles are controlled by parameters from Inputs data.

ctimes = s.iterate_state_to_convergence(cfc, m, grs, noplot=True)

After the convergence is reached or the maximum iteration number is reached, the method would return an instance of ComputingTime which contains time spent on each step of iteration. The user may plot or print the time data for performance review. Currently, there is no methods for post-processing results. The post-processing methods would be provide in the near future.

Performance

Current version uses Numba JIT (just-in-time) compiler to accelerate Python code. Thus, for the first run of ISOLVER, it will take a few extra seconds to compile the math utilities library. This part code will be ported into Cython in the following update so that it could call pre-compiled library directly without delays.

To get best compatibility, the multi-processing is disabled to get the code running with module “anaconda2” on PPPL cluster. Also the “fast math” option for Numba JIT is disabled for public interpreter “anaconda2”, since the “fast math” option uses Intel SVML to accelerate computing, which would lead to small errors for non-Intel platform.

Without plotting, the original IDL program would take 64 iterations to converge and cost ~80 seconds on AMD CPU or ~43 second on Intel CPU. The ported Python version would takes ~43 second on AMD CPU or ~26 second on Intel CPU. A better performance could be achieved with multi-processing and Intel SVML.

Documentation for the Code

LRDFIT

class lib.lrdfit.AuxCurrDr(ipwant=1.0, jip=1.0, jr0=None, jw=None)

Current Drive Parameters for Auxillary Plasma Heating

Parameters:

• ipwant (float, optional) – plasma current want (in Amperes)
• jip (float, optional) – current density
• jr0 (float, optional) –
• jw (float, optional) –

ip

float – plasma current (in Amperes)

r0

float

w

float

class lib.lrdfit.BetaCtrl(betamode=None, betascl=None, betatwant=None, betanwant=None)

Control parameters for the equilibrium pressure

Parameters:

• betamode (int, optional) –

  Beta mode:

  0 ==> beta-t = scaled initial equilibrium beta-t;

  1 ==> beta-t = betatwant;

  2 ==> beta-N = betanwant;

  3 ==> beta-N = from kinetic p profiles;

• betascl (float, optional) – scale factor for initial equilibrium beta
• betatwant (float, optional) – vacuum toroidal beta (in percent)
• betanwant (float, optional) – normalized beta (%mT / MA)

betamode

int – Beta mode:

0 ==> beta-t = scaled initial equilibrium beta-t;

1 ==> beta-t = betatwant;

2 ==> beta-N = betanwant;

3 ==> beta-N = from kinetic p profiles;

betascl

float – scale factor for initial equilibrium beta

betatwant

float – vacuum toroidal beta (in percent)

betanwant

float – normalized beta (%mT / MA)

class lib.lrdfit.BoundaryFitIcoil

Coil current data which has best fit to the reference surface

psit

np.ndarray – total poloidal flux

psic

np.ndarray – poloidal flux from all coils

icoil

np.ndarray – coil current

xaxis

np.ndarray – X coordinate of best-guess axis

zaxis

np.ndarray – Z coordinate of best-guess axis

psitaxis

np.ndarray – totol poloidal flux at axis

xref

np.ndarray – X coordinate of reference point on the boundary

zref

np.ndarray – Z coordinate of reference point on the boundary

psiref

np.ndarray – total poloidal flux at reference surface

rvec

np.ndarray – R coordinate of grid

zvec

np.ndarray – Z coordinate of grid

psipsgn

float – The sign of poloidal flux from plasma current

change_bfi_dimensions(nrnew, nznew)

Change dimensions of coil current data

Interpolate new data of new dimensions from existing coil current data, which has best fit to the reference surface.

Parameters:

• nrnew (int) – new r dimension
• nznew (int) – new z dimension

compute(xbref, zbref, psip, grs, wpsi, wrbp, icoil, wicoil)

Compute coil currents fitted to plasma boundary

Given the reference surface defined by (xbref,zbref) and the poloidal flux from the plasma (psip), compute the coil currents which best-fit a flux surface to the shape of the reference surface.

Parameters:

• xbref (np.ndarray) – X coordinates of the reference surface
• zbref (np.ndarray) – Z coordinates of the reference surface
• psip (np.ndarray) – poloidal flux from the plasma
• grs (GridStructure) – Grid data structure
• wpsi (np.ndarray) – weight of poloidal flux
• wrbp (np.ndarray) – weight of gradient of poloidal flux
• icoil (np.ndarray) – coil current
• wicoil (np.ndarray) – weight of coil current

class lib.lrdfit.CoilCurrParam(icoil=None, wicoil=None)

Coil current parameter structure

Parameters:

• icoil (np.ndarray, optional) – control coil current
• wicoil (np.ndarray, optional) – weighting factor for coil current

icoil

np.ndarray – control coil current

wicoil

np.ndarray – weighting factor for coil current

class lib.lrdfit.ComputingTime(tsize, t0=None)

Computing time used for each step

This class of data is structured to record computing time of each step in iterated convergence of LRDFIT data. Each attribute contains an array of time duration to compute that step for each iteration.

This class is used for benchmarking the performance of computing.

Parameters:

• tsize (int) – maximum number of iterations
• t0 (float, optional) – initial time

t1

np.ndarray – time interval to compute 1st step of each iteration

t2

np.ndarray – time interval to compute 2nd step of each iteration

t3

np.ndarray – time interval to compute 3rd step of each iteration

t4

np.ndarray – time interval to compute 4th step of each iteration

t5

np.ndarray – time interval to compute 5th step of each iteration

t6

np.ndarray – time interval to compute 6th step of each iteration

t7

np.ndarray – time interval to compute 7th step of each iteration

t8

np.ndarray – time interval to compute 8th step of each iteration

t9

np.ndarray – time interval to compute 9th step of each iteration

t10

np.ndarray – time interval to compute 10th step of each iteration

mean_time(tname)

Get the average computing time of any step of computing

Parameters:tname (str) – attribute name Returns:mean time Return type:float

print_time()

print mean execution time for each step

remove_empty_values(iters)

Remove unused empty time elements for each step

Parameters:iters (int) – number of total iterations

class lib.lrdfit.CurrCtrl(ipmode=None, ipscl=None, ipwant=None)

Control parameters for the equilibrium current

Parameters:

• ipmode (int, optional) –

  Plasma current mode:

  0 ==> Ip = scaled initial equilibrium Ip;

  1 ==> Ip = Ipwant;

  2 ==> Ip = from kinetic J profiles;

• ipscl (float, optional) – Scale factor for initial equilibrium Ip
• ipwant (float, optional) – Plasma current want

ipmode

int – Plasma current mode:

0 ==> Ip = scaled initial equilibrium Ip;

1 ==> Ip = Ipwant;

2 ==> Ip = from kinetic J profiles;

ipscl

float – Scale factor for initial equilibrium Ip

ipwant

float – Plasma current want

class lib.lrdfit.CurrPrflCtrl(ppsrc=None, jpsrc=None, ppmode=None, jpmode=None)

Plasma current profiles control parameter

Parameters:

• ppsrc (int, optional) –

  Control parameters for the equilibrium kinetic profile source:

  0 ==> p profile from initial equilibrium;

  1 ==> p profile from input parameterization;

• jpsrc (int, optional) –

  J profile source:

  0 ==> J profile from initial equilibrium;

  1 ==> J profile from input parameterization;

• ppmode (int, optional) –

  p profile mode:

  0 ==> Scale entire p to satisfy betamode;

• jpmode (int, optional) –

  J profile mode:

  0 ==> Scale entire J to satisfy ipmode;

  1 ==> Scale Vloop to satisfy ipmode when jpsrc > 0;

ppsrc

int – Control parameters for the equilibrium kinetic profile source:

0 ==> p profile from initial equilibrium;

1 ==> p profile from input parameterization;

jpsrc

int – J profile source:

0 ==> J profile from initial equilibrium;

1 ==> J profile from input parameterization;

ppmode

int – p profile mode:

0 ==> Scale entire p to satisfy betamode;

jpmode

int – J profile mode:

0 ==> Scale entire J to satisfy ipmode;

1 ==> Scale Vloop to satisfy ipmode when jpsrc > 0;

class lib.lrdfit.EqsPars(eqsrc=None, ipsign=None, btsign=None, isym=None, ibt0=None, r_bt0=None, bt0=None)

The parameters for equilibrium data of vacuum toroidal field.

Parameters:

• eqsrc (str, optional) – the source of equilibrium data
• ipsign (float, optional) – the sign of plasma current data
• btsign (float, optional) – the sign of toroidal field
• isym (bool, optional) – if force plasma flux to be up/down symmetric
• ibt0 (bool, optional) – if use initial toroidal field
• r_bt0 (float, optional) – reference major radius for vacuum toroidal field (m)
• bt0 (float, optional) – Toroidal field at r_bt0 (Tesla)

eqsrc

str – the source of equilibrium data

ipsign

float – the sign of plasma current data

btsign

float – the sign of toroidal field

isym

bool – if force plasma flux to be up/down symmetric

ibt0

bool – if use initial toroidal field

r_bt0

float – reference major radius for vacuum toroidal field (m)

bt0

float – Toroidal field at r_bt0 (Tesla)

Examples

>>> eqspars = EqsPars(eqsrc='efitmds',
...                   ipsign=1.0,
...                   btsign=-1.0,
...                   isym=True,
...                   ibt0=True,
...                   r_bt0=0.94,
...                   bt0=-1.0)

class lib.lrdfit.EquilibriumCtrl(psin_ieq=None, p_ieq=None, jb_ieq=None, ip_ieq=None, betat_ieq=None, betat_want=None, betan_want=None, ip_want=None)

Equilibrium control data

This is initial equilibrium profile and scalar data used to generate kinect profiles

Parameters:

• psin_ieq (np.ndarray, optional) – normalized initial equilibrium poloidal flux
• p_ieq (np.ndarray, optional) – initial equilibrium poloidal flux
• jb_ieq (np.ndarray, optional) – initial equilibrium <J.B>
• ip_ieq (np.ndarray, optional) – initial equilibrium plasma current
• betat_ieq (np.ndarray, optional) – initial equilibrium beta, which is the ratio of the plasma pressure to the magnetic pressure
• betat_want (np.ndarray, optional) – wanted beta
• betan_want (np.ndarray, optional) – wanted normalized beta, which is the normalized ratio of the plasma pressure to the magnetic pressure
• ip_want (float, optional) – wanted plasma current

psin_ieq

np.ndarray – normalized initial equilibrium poloidal flux

p_ieq

np.ndarray – initial equilibrium poloidal flux

jb_ieq

np.ndarray – initial equilibrium <J.B>

ip_ieq

np.ndarray – initial equilibrium plasma current

betat_ieq

np.ndarray – initial equilibrium beta, which is the ratio of the plasma pressure to the magnetic pressure

betat_want

np.ndarray – wanted beta

betan_want

np.ndarray – wanted normalized beta, which is the normalized ratio of the plasma pressure to the magnetic pressure

ip_want

float – wanted plasma current

class lib.lrdfit.ExternalLibrarySignature(im, name, bits)

Signature for external shared library

Parameters:

• im (str) – the file path of shared object
• name (str) – the name of external library
• bits (int) – the type of processor 64 bits or 32 bits

class lib.lrdfit.FieldGrid

Field grid data structure

rmin

np.ndarray – the minimum major radius value on the grid

rmax

np.ndarray – the maximum major radius value on the grid

zmin

np.ndarray – the minimum vertical position on the grid

zmax

np.ndarray – the maximum vertical position on the grid

nr

int – the number of radial nodes for the low-resolution grid

nz

int – the number of vertical nodes for the low-resolution grid

nrh

int – the number of radial nodes for the high-resolution grid

nzh

int – the number of vertical nodes for the high-resolution grid

lims

Limiter – limiter coordinate data

class lib.lrdfit.FluxCoordParam(radexp0=None, radexp1=None, fluxfrac=None, nradius=101, ntheta=101, rscale=1.0, dconv=0.001, miters=1)

Parameters to solve flux coordinates

radexp0

float, optional – psi scales as rminor**radexp0 near axis

radexp1

float, optional – psi scales as rminor**radexp1 near edge

fluxfrac

float, optional – Fraction of total flux spanned by surfaces

nradius

int, optional – Number of flux surfaces

ntheta

int, optional – Number of theta points on surface

rscale

float, optional – Extrapolate past LCFS this far when interpolating

dconv

float, optional – Flux surface convergence tolerance

miters

int, optional – Maximum number of iterations to solve flux coordinates

Parameters:

• radexp0 (float) – psi scales as rminor**radexp0 near axis
• radexp1 (float) – psi scales as rminor**radexp1 near edge
• fluxfrac (float) – Fraction of total flux spanned by surfaces
• nradius (int) – Number of flux surfaces
• ntheta (int) – Number of theta points on surface
• rscale (float) – Extrapolate past LCFS this far when interpolating
• dconv (float) – Flux surface convergence tolerance
• miters (int) – Maximum number of iterations to solve flux coordinates

class lib.lrdfit.FluxSurfParam(nbndry=None, nbfit1=None, nbfit2=None)

Last closed flux surface parameter structure

Parameters:

• nbndry (int, optional) – Number of theta points on last closed flux surface (LCFS)
• nbfit1 (int, optional) – Number of radial ray points in finding plasma LCFS
• nbfit2 (int, optional) – Number of radial ray points in refining plasma LCFS

nbndry

int – Number of theta points on last closed flux surface (LCFS)

nbfit1

int – Number of radial ray points in finding plasma LCFS

nbfit2

int – Number of radial ray points in refining plasma LCFS

class lib.lrdfit.FreeBoundary

Free boundary data

xbnd

np.ndarray – X coordinate of free boundary

zbnd

np.ndarray – Z coordinate of free boundary

xaxis

np.ndarray – X coordinate of the axis

zaxis

np.ndarray – Z coordinate of the axis

psibnd

np.ndarray – total poloidal flux at free boundary

psiaxis

np.ndarray – total poloidal flux at the axis

thetabnd

np.ndarray – angle array with uniform step length

dpsidrb

np.ndarray – derivative of total poloidal flux with respect to R coordinate at free boundary

dpsidzb

np.ndarray – derivative of total poloidal flux with respect to Z coordinate at free boundary

find(bfi, grs, nbdry=101, nfit1=101, nfit2=11, xpt=None)

Find the new free boundary based on coil current data

Given the result of the boundary fit, find the new free boundary. The reference flux from “bfi” is used for the boundary flux.

Parameters:

• bfi (BoundaryFitIcoil) – coil currents fitted to plasma boundary
• grs (GridStructure) – Grid data structure
• nbdry (int, optional) – number of points on boundary
• nfit1 (int, optional) – number of points on rays from axis for finding boundary
• nfit2 (int, optional) – once boundary is bracketed, use this many subdivisions.
• xpt (PoloidalFluxProfile, optional) – x-point profile for free boundary

class lib.lrdfit.GEQPars(geqshot=None, geqtime=None, geqindex=None, geqexpt=None, geqname=None)

G-EQDSK parameters

A standard format for storing tokamak equilibrium data is G-EQDSK [1] which contains the poloidal flux in (R,Z) and 1D profiles of pressure, f = R * B(phi), safety factor q, and other quantities related to the solution. This class contains parameters for G-EQDSK or EFIT-MDS input data

Parameters:

• geqshot (int, optional) – the shot number from G-EQDSK file
• geqtime (list of floats, optional) – the time of shot from G-EQDSK file
• geqindex (int, optional) – the index from G-EQDSK file
• geqexpt (str, optional) – the name of experiment from G-EQDSK file
• geqname (str, optional) – the name and path of G-EQDSK file

geqshot

int – the shot number from G-EQDSK file

geqtime

list of floats, optional – the time of shot from G-EQDSK file

geqindex

int – the index from G-EQDSK file

geqexpt

str – the name of experiment from G-EQDSK file

geqname

str – the name and path of G-EQDSK file

Examples

>>> geqpars = GEQPars(geqshot=116313,
...                   geqtime=[0.860],
...                   geqindex=-6,
...                   geqexpt='nstx',
...                   geqname='geqdsk/g121241.00270.NgI103QG')

References

[1]	G-EQDSK file format documentation.

generate_jsolver_name(outextn)

generate name string for J-SOLVER

Parameters:outextn (str) – output extension Returns:the output name of jsolver Return type:str

class lib.lrdfit.GridStructure

Grid data structure

Contains grid structure data Refer to IDL subroutine: generate_grid_structure() in lrdfit_routines.pro

nr

int – number of R coordinates

nz

int – number of Z coordinates

nr1

int – number of R coordinates - 1

nz1

int – number of Z coordinates - 1

nrnz

int – number of R coordinates * number of Z coordinates

nr1nz1

int – (number of R coordinates - 1) * (number of Z coordinates - 1)

ncoils

int – number of coils

nchars

int – number of characters for coil names

coil_names

list of str – a list of coil names

coil_width

np.ndarray – coil widths

coil_height

np.ndarray – coil heights

coil_nturns

np.ndarray – number of coil turns

coil_caf

np.ndarray – coil caf

rvec

np.ndarray – R coordinate on generated uniform linear space

zvec

np.ndarray – Z coordinate on generated uniform linear space

rv

np.ndarray – R coordinate on mesh grid vertices as a 1D list

zv

np.ndarray – Z coordinate on mesh grid vertices as a 1D list

rva

np.ndarray – R coordinate on mesh grid as 2D array

zva

np.ndarray – Z coordinate on mesh grid as 2D array

rc

np.ndarray – R coordinate of mesh center point as a 1D list

zc

np.ndarray – Z coordinate of mesh center point as a 1D list

rca

np.ndarray – R coordinate of mesh center point as a 2D array

zca

np.ndarray – Z coordinate of mesh center point as a 2D array

rlim

np.ndarray

zlim

np.ndarray

psic

np.ndarray – poloidal flux of coil on the grid for each coil group

change_grs_dimensions(nrnew, nznew)

change the dimensions of GridStructure function data

Parameters:

• nrnew (int) – new dimension of r
• nznew (int) – new dimension of z

generate_grid_structure(grddat, hr=False)

Generate grid data structure from limiter data and field grid data

From the original IDL “function generate_grid_structure, limdat, grddat, hr=hr” In the new function, “limdat” is replaced by an attribute of FieldGrid().lims

Parameters:

• grddat (FieldGrid) – field grid data
• hr (bool, optional) – if use high resolution grid

get_greens_data(filename)

Get grid data from file for Green’s function

Get Grid data structure from existing data file. Those data would be used to solve differential equations using Green’s method.

Parameters:filename (str) – file name and path

class lib.lrdfit.Inputs

ISOLVER input data and parameters

input_count

int – counter for how many times the input has been changed

eqspars

EqsPars – Applied vacuum TF parameters

jsopars

JsoPars – Shot and time parameters for OUTPUT GEQDSK file

trpars

TRANSPEqm – Parameters for TRANSP input data

geqpars

GEQPars – Parameters for GEQDSK or EFIT-MDS input data

betactls

BetaCtrl – Control parameters for the equilibrium pressure

ipctls

CurrCtrl – Control parameters for the equilibrium current

kpctls

CurrPrflCtrl – Plasma current profiles control parameters

lrdctls

LRDFITCtrl – Parameters for writing ISOLVER output files to LRDFIT tree

itconvpars

IterConvParam – iteration convergence parameters

lcfspars

FluxSurfParam – Free boundary equilibrium solution parameters

limpars

Limiter – Limiter surface parameters

fcspars

FluxCoordParam – Fixed boundary flux-coordinate equilibrium solution parameters

shapepars

ShapeParam – All shape parameters for shape control

spcpars

StrikePointCtrl – Strike point control parameters

icoilpars

CoilCurrParam – Coil current control parameters

kprofpars

KineticProfileParam – kinetic profile parameters

load_json(filename)

Read JSON file into current input instance

Parameters:filename (str) – file name and path to read JSON file

save_json(filename)

Save current input object as JSON file

Parameters:filename (str) – file name and path to write JSON file

class lib.lrdfit.IterConvParam(nrgrid=None, nzgrid=None, ngdoub=None, nitmax=None, newpsifrc=None, convtol=0.001, convctl=None)

Iteration convergence parameters

Parameters:

• nrgrid (int, optional) – Initial number of major radial (R) grid points of psi(R,Z) solution grid
• nzgrid (int, optional) – Initial number of vertical (Z) grid points of psi(R,Z) solution grid
• ngdoub (int, optional) – Number of grid resolution (R & Z) doublings during convergence cycle
• nitmax (int, optional) – Maximum number of iterations allowed
• newpsifrc (float, optional) – Fraction of new plasma poloidal flux used during iteration
• convtol (float, optional) – Convergence tolerance
• convctl (np.ndarray, optional) – Control parameters for convergence tolerance of each profile

nrgrid

int – Initial number of major radial (R) grid points of psi(R,Z) solution grid

nzgrid

int – Initial number of vertical (Z) grid points of psi(R,Z) solution grid

ngdoub

int – Number of grid resolution (R & Z) doublings during convergence cycle

nitmax

int – Maximum number of iterations allowed

newpsifrc

float – Fraction of new plasma poloidal flux used during iteration

qprofctol

float – q profile convergence tolerance

fprofctol

float – F profile convergence tolerance

pprofctol

float – p profile convergence tolerance

icoilctol

float – Coil current convergence tolerance

frbndctol

float – Free boundary spatial convergence tolerance (meters)

prssrctol

float – Plasma pressure (beta-N) convergence tolerance

class lib.lrdfit.JPHIS

Toroidal current density data

Jtorv

np.ndarray – Current density on full grid at torodial vertices

Itorc

np.ndarray – Current on sub- grid at torodial centers (Amperes)

ip

float – plasma current

dr

float – the size of small change in R direction

dz

float – the size of small change in Z direction

static call_external_psip(_a, _b, _c, _d, _f, _m, _n)

Wrapper function to call the external FORTRAN library

Solve for the poloidal flux psi = R * Aphi from the plasma current given the plasma toroidal current density Jphi on the solution grid

Parameters:

• lsod (ExternalLibrarySignature) – shared object path
• _a (np.ndarray) – Rmin on rectangular, uniformly spaced, R,Z grid
• _b (np.ndarray) – Rmax on rectangular, uniformly spaced, R,Z grid
• _c (np.ndarray) – Zmin on rectangular, uniformly spaced, R,Z grid
• _d (np.ndarray) – Zmax on rectangular, uniformly spaced, R,Z grid
• _f (np.ndarray) – Jphi in MKS units, m1 x n1 2D array
• _m (np.int) – number of major radial intervals on grid
• _n (np.int) – number of vertical intervals on grid

Returns:

if runs OK, return would be 0

Return type:

int

change_jphis_dimensions(nrnew, nznew)

change the dimensions of JPHIS object

Parameters:

• nrnew (int) – new dimension of r
• nznew (int) – new dimension of z

jphis_eqs(eqs, grs, ipsign=1.0, ipwant=None, fast=None)

Derive jphis data from EQS data

Parameters:

• eqs (eqsc.Equilibrium) – equilibrium solution
• grs (GridStructure) – grid data structure
• ipsign (float, optional) – the sign of plasma current
• ipwant (float, optional) – the plasma current wanted
• fast (bool, optional) – fast method to get approximate current density

jphis_jphi(jphi, grs, ipsign=None)

Generate the Jphi-structure from Jphi on the primary grid and grs

Parameters:

• jphi (np.ndarray) – Jphi on the primary grid
• grs (GridStructure) – GridStructure function data
• ipsign (float, optional) – sign of Ip

psip_jphis(grs)

Compute the poloidal flux on grid from plasma current density on grid

This is a wrapper for internal method to call external library

Parameters:grs (GridStructure) – Grid data structure

class lib.lrdfit.JsoPars(jsonamein=None, geqpars=None, outextn=None)

JSOLVER parameters

JSOLVER is a fixed-boundary, toroidal MHD equilibrium code, which solves the Grad-Shafranov equation for the poloidal flux using a finite difference scheme. JSOLVER applies an iterative metric method to simultaneously solve for the plasma equilibrium and the magnetic in flux in equal-arc coordinates. This class contains shot and time parameters for OUTPUT GEQDSK file from JSOLVER.

Parameters:

• jsonamein (str, optional) – the input file name from JSOLVER
• geqpars (GEQPars, optional) – the parameters of G-EQDSK data
• outextn (str, optional) – the filename extension for G-EQDSK output file

jsonamein

str – the input file name from JSOLVER

jsonameout

str – the output file name from JSOLVER

class lib.lrdfit.KineticCurrentProfile

Kinetic current profile

This class inherits from neoc.BootstrapCurrent class with extra attributes

psin

np.ndarray – normalized poloidal flux

jkin

np.ndarray – kinetic current density

ip_jkin

np.ndarray – kinetic plasma current

ip_ktot

np.ndarray – total kinetic plasma current

jauxb

np.ndarray – beam current density generated by auxillary current drive

jauxc

np.ndarray – core current density generated by auxillary current drive

jauxm

np.ndarray – mid-radius current density generated by auxillary current drive

jaux

np.ndarray – total current density generated by auxillary current drive

ip_auxb

np.ndarray – beam current generated by auxillary current drive

ip_auxc

np.ndarray – core current generated by auxillary current drive

ip_auxm

np.ndarray – mid-radius current generated by auxillary current drive

ip_aux

np.ndarray – total current generated by auxillary current drive

derive_from_state(s, c, grs)

Derive kinetic current profile from an instance of class State and other data

Parameters:

• s (State) – tokamak state data
• c (FluxCoord) – flux coordinate data
• grs (GridStructure) – Grid data structure

class lib.lrdfit.KineticParam(form=None, kp0=1.0, kp1=1.0, w1=None, a1=None, b1=None, a2=None, b2=None, etfc=None, sc=None)

Kinetic parameter for plasma

Parameters:

• form (int, optional) –

  Profile functional form index:

  form == 1 ==> # extended tanh-function;

  form == 2 ==> # transport code - sqrt(psitor);

• etfc (np.ndarray, optional) – Extended tanh-function
• sc (float, optional) – scale factor
• kp0 (float, optional) – Kinetic parameter 0
• kp1 (float, optional) – Kinetic parameter 1
• w1 (float, optional) – Weight parameter 1
• a1 (float, optional) – Profile parameter a1
• b1 (float, optional) – Profile parameter b1
• a2 (float, optional) – Profile parameter a2
• b2 (float, optional) – Profile parameter b2

form

int – Profile functional form index:

form == 1 ==> # extended tanh-function;

form == 2 ==> # transport code - sqrt(psitor);

etfc

np.ndarray – Extended tanh-function

kp0

float – Kinetic parameter 0

kp1

float – Kinetic parameter 1

w1

float – Weight parameter 1

a1

float – Profile parameter a1

b1

float – Profile parameter b1

a2

float – Profile parameter a2

b2

float – Profile parameter b2

class lib.lrdfit.KineticProfile(kev=None, cc=None)

Kinetic profile

Parameters:

• kev (bool, optional) – if True, use ‘keV’ as unit for electron temperature, otherwise, use ‘eV’
• cc (bool, optional) – if True, use ‘cm^-3’ as unit for density, otherwise, use ‘m^-3’

tscale

float – scale of electron temperature profile

nscale

float – scale of square root of normalized profile

vscale

float – scale of voltage profile

tuscale

float – scale of electron temperature unit

nuscale

float – scale of normalized profile unit

nunits

str – unit of normalized profile

tunits

str – unit of electron temperature

dunits

str – unit of density

punits

str – unit of pressure

junits

str – unit of current density

vunits

str – unit of voltage

rhotor

KineticProfileStruct – kinetic profile of square root of normalized psi_tor, psi_tor = t_flux / t_fluxbdy where t_flux = the toroidal flux enclosed by the surface, and t_fluxbdy is the total toroidal flux enclosed within the surface defining the plasma boundary.

rhovol

KineticProfileStruct – kinetic profile of square root of normalized volume

rhopol

KineticProfileStruct – kinetic profile of square root of normalized psi_pol, which is the poloidal flux enclosed within the surface.

psin

KineticProfileStruct – kinetic profile of normalized poloidal flux

Te

KineticProfileStruct – kinetic profile of electron temperature

Ti

KineticProfileStruct – kinetic profile of main ion temperature

Tz

KineticProfileStruct – kinetic profile of impurity ion temperature

Tf

KineticProfileStruct – kinetic profile of fast ion effective temperature

ne

KineticProfileStruct – kinetic profile of electron density

ni

KineticProfileStruct – kinetic profile of main ion density

nz

KineticProfileStruct – kinetic profile of impurity ion density

nf

KineticProfileStruct – kinetic profile of fast ion density

Za

KineticProfileStruct – kinetic profile of average charge

Ze

KineticProfileStruct – kinetic profile of effective charge

pe

KineticProfileStruct – kinetic profile of electron pressure

pi

KineticProfileStruct – kinetic profile of main ion pressure

pz

KineticProfileStruct – kinetic profile of impurity ion pressure

pf

KineticProfileStruct – kinetic profile of fast ion pressure

pt

KineticProfileStruct – kinetic profile of thermal ion pressure

p

KineticProfileStruct – kinetic profile of total pressure

Jb

KineticProfileStruct – kinetic profile of neutral beam injection current density

Jc

KineticProfileStruct – kinetic profile of core-localized auxiliary current density

Jm

KineticProfileStruct – kinetic profile of mid-radius auxiliary current density

Vl

KineticProfileStruct – kinetic profile of loop voltage

derive_from_kprofpars(inputobj, psin=None, nscale=1.0, tscale=1.0, vscale=1.0, m=None, tcd=None)

Derive kinectic profiles from kinetic profile parameters

Parameters:

• inputobj (Inputs) – input data
• psin (np.ndarray, optional) – the poloidal flux enclosed within the surface
• nscale (float, optional) – scale factor of n
• tscale (float, optional) – scale factor of T
• vscale (float, optional) – scale factor of V
• m (Momentum, optional) – state moment data
• tcd (KineticProfile, optional) – kinetic profiles from transport code

derive_from_state(s, m)

Derive kinetic profiles from current state

Parameters:

• s (State) – an instance of State class
• m (Momentum) – state of moments

class lib.lrdfit.KineticProfileParam(nscale=None, tscale=None, vscale=None, jb_form=None, newjbf=None, minjbr=None, vlmin=None, vlmax=None)

Kinetic profile parameter structure

Parameters:

• nscale (float, optional) – Kinetic profile scales of n
• tscale (float, optional) – Kinetic profile scales of T
• vscale (float, optional) – Kinetic profile scales of V
• jb_form (int, optional) – Jaux - NBI
• newjbf (float, optional) – Mix new <J.B> solution into old with this fraction
• minjbr (float, optional) – threshold to avoid current hole
• vlmin (float, optional) – Loop voltage limits - minimum
• vlmax (float, optional) – Loop voltage limits - maximum

nscale

float – Kinetic profile scales of n

tscale

float – Kinetic profile scales of T

vscale

float – Kinetic profile scales of V

jb_form

int – Jaux - NBI

jnbi

AuxCurrDr – J - NBI current

jcen

AuxCurrDr – Jaux - plasma core

jmid

AuxCurrDr – Jaux - mid radius

newjbf

float – Mix new <J.B> solution into old with this fraction

minjbr

float – threshold to avoid current hole

vlmin

float – Loop voltage limits - minimum

vlmax

float – Loop voltage limits - maximum

Z

Species – species charges

m

Species – species masses

Te

KineticParam – Electron temperature kinetic parameters

ne

KineticParam – Electron density kinetic parameters

Ti

KineticParam – Thermal ion temperature kinetic parameters

Ze

KineticParam – Effective ion charge kinetic parameters

nf

KineticParam – Fast ion density kinetic parameters

Tf

KineticParam – Fast ion temperature kinetic parameters

Vl

KineticParam – Loop voltage kinetic parameters

generate_profile(psink, profname, rhotor, tcd=None)

generate kinetic profile

Generate kinetic profile from preset parameters, extended tanh-function or kinetic profile from transport code.

Parameters:

• psink (np.ndarray) – The kinetic poloidal flux enclosed within the surface
• profname (str) – Profile name
• rhotor (np.ndarray) – Square root of normailized psi_tor
• tcd (KineticProfile, optional) – Kinetic profiles from transport code

Returns:

Generated kinectic profile

Return type:

np.ndarray

Examples

>>> ndrho = 201
>>> drho = np.arange(ndrho, dtype=float) / (ndrho - 1.0)
>>> psink = drho**2
>>> rhotor = np.zeros(ndrho)
>>> tcd = KineticProfile()
>>> profile = self.generate_profile(psink, 'Te', rhotor, tcd=tcd)

class lib.lrdfit.KineticProfileStruct(desc=None, units=None, profile=None, data=None)

Structured data of kinetic profiles

Parameters:

• desc (str, optional) – description of the kinetic profile
• units (str, optional) – unit of the kinetic data
• profile (np.ndarray, optional) – profile data
• data (np.ndarray, optional) – kinetic data

description

str – description of the kinetic profile

units

str – unit of the kinetic data

profile

np.ndarray – profile data

data

np.ndarray – kinetic data

class lib.lrdfit.LRDFITCtrl(shot=None, index=None, wrtpath=None, time0=None, dt=None, nt=None)

Parameters for writing ISOLVER output files to LRDFIT (MDSplus) tree

Parameters:

• shot (int, optional) – LRDFIT shot number
• index (float, optional) – LRDFIT index
• wrtpath (float, optional) – LRDFIT wrt path
• time0 (int, optional) – LRDFIT time 0
• dt (int, optional) – LRDFIT dt
• nt (str, optional) – LRDFIT nt

shot

int – LRDFIT shot number

eindex

float – LRDFIT index

wrtpath

float – LRDFIT wrt path

time0

int – LRDFIT time 0

dt

int – LRDFIT dt

nt

str – LRDFIT nt

Examples

>>> lrdctls = LRDFITCtrl(200160, 0.0, 0.1, 11, 1, nt='./fits')

class lib.lrdfit.Limiter(nlimh=None)

Limiter contour data

The Limiter data contains R,Z coordinates, which should form a closed curve in the poloidal plane.

A limiter is a material surface within the tokamak vessel which defines the edge of the plasma and thus avoids contact between the plasma and the vessel. Alternative to using a divertor to define the edge.

Parameters:nlimh (int, optional) – Number of (R,Z) points on high-resolution limiter boundary

rlim

np.ndarray – limiter R coordinates form a closed curve in the poloidal plane

zlim

np.ndarray – limiter Z coordinates form a closed curve in the poloidal plane

rlimh

np.ndarray – R coordinates on the high-resolution limiter surface

zlimh

np.ndarray – Z coordinates on the high-resolution limiter surface

ivilim

np.ndarray – indices of grid vertices enclosed by limiter boundary

icilim

np.ndarray – indices of grid centers enclosed by limiter boundary

ncilim

int – number of grid centers enclosed by limiter boundary

nvertex

int – number of vertices required for a zone to be inside the limiter

nlimh

int – Number of (R,Z) points on high-resolution limiter boundary

nsubdiv

int – the number of subdivisions allowed in refining the size of the regions for plasma elements near the limiter boundary

limiter_data_structure(grs, inputobj, use_nvertex=None)

Derive limiter data from Input data and GridStructure data

Parameters:

• grs (GridStructure) – Grid data structure
• inputobj (Input) – an instance of class Inputs
• use_nvertex (bool, optional) – use input parameters for how many number of vertices must be enclosed to count region as enclosed

class lib.lrdfit.OutParam(geqextn=None)

Output parameters for GEQDSK file

Parameters:geqextn (str, optional) – equilibrium output extension name

geqextn

str – equilibrium output extension name

class lib.lrdfit.PlasmaBoundary

Plasma boundary data including x-points and strike points

X-point is a point at which core plasma and scrape-off layer plasma meet on the separatrix.

Strike points are points at where the particles escaping from the plasma and following the field lines outside the seperatrix, collide with the walls.

newfb

FreeBoundary – new free boundary data

newbfi

BoundaryFitIcoil – new coil currents fitted to plasma boundary

usexptb

bool – determine if use x-point boundary to compute strike-point coordinates

nxpt

int – number of x-points

xptb

PoloidalFluxProfile – poloidal flux profile of x-point at boundary

xptv

PoloidalFluxProfile – poloidal flux profile of x-point array

xpt1

PoloidalFluxProfile – poloidal flux profile of lower x-point

xpt2

PoloidalFluxProfile – poloidal flux profile of upper x-point

spc

StrikePointCoords – strike point data

drsep

float – radial distance at the outer midplane between the flux surfaces connected to the upper and lower x-points

compute(newbfi, grs, lims, nbdry=None, nfit1=None, nfit2=None, quiet=None, fedge=None, debug=False)

Compute plasma boundary data

Parameters:

• newbfi (BoundaryFitIcoil) – coil currents fitted to plasma boundary
• grs (GridStructure) – Grid data structure
• lims (Limiter) – Limiter contour data
• nbdry (int, optional) – dimension of boundary data
• nfit1 (int, optional) – number of points on rays from axis for finding boundary
• nfit2 (int, optional) – once boundary is bracketed, use this many subdivisions.
• quiet (bool, optional) – option to turn on/off debug message
• fedge (np.ndarray, optional) – force at plasma edge (?)
• debug (bool, optional) – option to turn on/off debug mode

compute_state_drsep(quiet=None)

Compute the drsep of plasma boundary

Parameters:quiet (bool, optional) – option to turn on/off debug message

class lib.lrdfit.PoloidalFluxProfile(r=None, z=None, psi=None, i=None)

Poloidal profile

This includes the coordiantes and poloidal flux.

Parameters:

• r (np.ndarray, optional) – R coordinate
• z (np.ndarray, optional) – Z coordinate
• psi (np.ndarray, optional) – poloidal flux
• i (int, optional) – index of points

r

np.ndarray – R coordinate

z

np.ndarray – Z coordinate

psi

np.ndarray – poloidal flux

i

int – index of point

class lib.lrdfit.ReferenceShape(x, z)

Boundary shape data

Parameters:

• x (np.ndarray) – boundary data in x coordinate
• z (np.ndarray) – boundary data in z coordinate

interpolate(nt)

Interpolate the reference shape

This is part of unfinished option which allows user to use mix of experimental and user-specified boundary shapes

Parameters:nt (int) – dimension of theta

class lib.lrdfit.ShapeControlReference

Reference data for plasma shape control

x

np.ndarray – X coordinates of the reference surface

z

np.ndarray – Z coordinates of the reference surface

wpsi

np.ndarray – weight of poloidal flux

wrbp

np.ndarray – weight of gradient of poloidal flux

theta

np.ndarray – poloidal angle

tnorm

np.ndarray – normalized poloidal angle 0.0 <= tnorm < 1.0

define_shape_control_reference(inputobj, shape=None)

Define the reference surface for shape control

Parameters:

• inputobj (Inputs) – input data for ISOLVER
• shape (ReferenceShape, optional) – allow user to use user specified boundary shape

class lib.lrdfit.ShapeCtrlParam(upper=None, lower=None)

Shape control parameters

Parameters:

• upper (optional) – Parameter for upper bound
• lower – Parameter for lower bound

upper

optional – Parameter for upper bound

lower

Parameter for lower bound

class lib.lrdfit.ShapeParam(rshift=None, zshift=None, scminorrad=None, minorrad=None, r0=None, z0=None, nshape=None)

Shape parameters

Parameters:

• rshift (float, optional) – Shift of the JSOLVER boundary in R axis
• zshift (float, optional) – Shift of the JSOLVER boundary in Z axis
• scminorrad (float, optional) – scale of minor radius
• minorrad (float, optional) – plasma minor radius
• r0 (float, optional) – major radius of plasma geometric center
• z0 (float, optional) – vertical position of plasma geometric center
• nshape (int, optional) – number of shape control points

rshift

float – Shift of the JSOLVER boundary in R axis

zshift

float – Shift of the JSOLVER boundary in Z axis

newbnd

ShapeCtrlParam – Option for new boundary:

newbnd.upper == 0 ==> use supplied reference shape;

newbnd.upper == 1 ==> use parameterized JSOLVER boundary;

scminorrad

float – scale of minor radius

minorrad

float – plasma minor radius

r0

float – major radius of plasma geometric center

z0

float – vertical position of plasma geometric center

nshape

int – number of shape control points

k

ShapeCtrlParam – elongation

d

ShapeCtrlParam – triangularity

zo

ShapeCtrlParam – outer squareness

zi

ShapeCtrlParam – inner squareness

wo

ShapeCtrlParam – outer squareness weight during flux control

wi

ShapeCtrlParam – inner squareness weight during flux control

xpt

ShapeCtrlParam – x-point control parameter xpt.upper == 1 ==> generate x-point

wxpt

ShapeCtrlParam – x-point weighting relative to flux control

class lib.lrdfit.SharedObject

Information about JPHI shared library.

The FORTRAN library to solve Grad-Shafranov equations using Von Hagenov method. There are three functions from the shared library:

• psip_from_jphi
• psib_from_jphi
• aphi_from_jphi

psip_from_jphi

ExternalLibrarySignature – Solve for the poloidal flux psi = R * Aphi from the plasma current given the plasma toroidal current density Jphi on the solution grid

psib_from_jphi

ExternalLibrarySignature – Compute flux on boundary using Von Hagenow matrix and dU/dn array

aphi_from_jphi

ExternalLibrarySignature – This routine computes the toroidal component of the vector potential in cylindrical (R,phi,Z) coordinates given the toroidal component of the current density on the same grid. The routine uses HWSCYL, which

class lib.lrdfit.Species(p=None, f=None, i=None)

Species of charges or masses

Parameters:

• p (int, optional) – species of bulk plasma
• f (int, optional) – species of fast ion
• i (int, optional) – species of impurity

p

int – species of bulk plasma

f

int – species of fast ion

i

int – species of impurity

class lib.lrdfit.State(inputobj)

Current state of teh 2D axisymmetric circuit model of a tokamak

This class is the master class of LRDFIT to store the state of each iteration and includes methods to initialize and converge the model data.

Parameters:inputobj (Inputs) – input data to initialize state

inputobj

Inputs – input data to initialize state

jphis

JPHIS – Toroidal current density data

psip

np.ndarray – poloidal flux on grid generated from the plasma

bfi

BoundaryFitIcoil – coil current data which has best fit to the reference surface

pb

PlasmaBoundary – plasma boundary data including x-points and strike points

lims

Limiter – limiter contour data

eqctls

EquilibriumCtrl – Equilibrium control data

rprofs

StateProfile – reference state profle

fprofs

StateProfile – fitted state profile

tcdata

KineticProfile – kinetic profile from transport code data

kprofs

KineticProfile – kinetic profile data

jprofs

KineticCurrentProfile – kinetic current profile

rshape

ShapeControlReference – Reference data for plasma shape control

figure

plt.Figure – the handler for plotting figure and set format

__contour0

plt.Contour – handler to update a set of countours for poloidal flux

__contour1

plt.Contour – if use x-point boundary to compute strike-point coordinates, then, use this handler to update a set of contours

change_state_grid_dimensions(grs, nrnew, nznew)

Change the grid dimensions of current state data

Parameters:

• grs (GridStructure) – Grid data structure
• nrnew (new r dimension size) –
• nznew (new z dimension size) –

compute_state_convergence_error(ostate, omoments, nmoments, cfc)

Compute the convergence error between new state and old state

Parameters:

• ostate (State) – state data in previous iteration
• omoments (Momentum) – state moment data in previous iteration
• nmoments (Momentum) – state moment data in current iteration
• cfc (FluxCoord) – flux coordinate data in current iteration

compute_state_jphi(psin, r_psin)

Use the reference flux functions to compute Jphi at specified values of psin and R(psin) for psin less than 1

Parameters:

• psin (np.ndarray) – normalized poloidal flux
• r_psin (np.ndarray) – r coordinate of normalized poloidal flux

Returns:

toroidal current density as sum of three terms

Return type:

float

double_state_grid_dimensions(grs)

Double the dimensiosn of grid for an instance of class State

Parameters:grs (GridStructure) – Grid data structure

initialize_from_equilibrium(grs, cfc, fast=None, psifrac=None)

Initialize current state from equilibrium data

Parameters:

• grs (GridStructure) – Grid data structure
• cfc (FluxCoord) – flux coordinate data
• fast (bool, optional) – an option to use fast method
• psifrac (float, optional) – fraction of poloidal flux spanned by flux surfaces

iterate_state_to_convergence(fc, m, grs, noplot=None)

Iterate the state profiles to convergence

Parameters:

• fc (FluxCoord) – flux coordinate data
• m (Momentum) – state moment data
• grs (GridStructure) – Grid data structure
• noplot (bool, optional) – option to turn on/off the plotting method

Returns:

computing time for each step of iteration

Return type:

ComputingTime

plot_state_convergence(s, fc, grs, plotbcp=None)

update plot figures about convergence between previous state and new state

Parameters:

• s (State) – previous State before the iteration
• fc (FluxCoord) – flux coordinate data
• grs (GridStructure) – Grid data structure
• plotbcp (bool, optional) – option to control if plot bcp

plot_state_init(grs)

initialize plot for state

initialize empty canvas for plotting convergence between previous state and new state

Parameters:grs (GridStructure) – Grid data structure

rescale_state_profiles(moments, printm=True)

Rescale equilibrium profiles

Rescale the euilibrium profiles of p, p’’ and <J.B>/<Bt/R2> and print the moment data and current profiles.

Parameters:

• moments (Momentum) – state moment data to derive scale for equilibrium profiles
• printm (bool, optional) – option to print key values of moment data

rescale_toroidal_field()

rescale toroidal field

rescale the reference toroidal field profile

update_state_boundary(grs)

Update plasma boundary data in current state

Parameters:grs (GridStructure) – Grid data structure

update_state_jphis(grs)

Update toroidal current density data

Parameters:grs (GridStructure) – Grid data structure

update_state_kinetic_profiles(cfc, grs, m)

Update the state of kinetic profiles in current state

Parameters:

• cfc (FluxCoord) – flux coordinate data
• grs (GridStructure) – Grid data structure
• m (Momentum) – state moment data

update_state_profiles(cfc)

Update state profile

update the profile of tokamak state based on the flux coordinate data

Parameters:cfc (FluxCoord) – flux coordinate data

update_state_psip(grs, psi_passive=None)

Update poloidal flux generated from the plasma

Parameters:

• grs (GridStructure) – Grid data structure
• psi_passive (np.ndarray, optional) – externally supplied flux from passive structure

update_state_totalflux(grs, filepath=None)

Update total flux data

Parameters:

• grs (GridStructure) – Grid data structure
• filepath (str, optional) – if filepath is None, print coil current information if filepath is not None, write coil current information to local file

class lib.lrdfit.StateProfile(psi=None, psin=None, pp=None, ffp=None, p=None, f=None, jb=None, q=None, rhovn=None, epot=None, m2=None, m2p=None, bt2=None, bp2=None, b2=None, fedge=None)

State profile

The profile class is to construct state profile data in LRDFIT

Parameters:

• psi (np.ndarray, optional) – poloidal flux
• psin (np.ndarray, optional) – normalized poloidal flux
• ffp (np.ndarray, optional) – FF’ profile
• f (np.ndarray, optional) – F profile
• pp (np.ndarray, optional) – p and p’ profile
• p (np.ndarray, optional) – p profile
• jb (np.ndarray, optional) – <J.B> profile
• q (np.ndarray, optional) – q profile
• rhovn (np.ndarray, optional) – normalized rho
• epot (np.ndarray, optional) – electrostatic potential
• m2 (np.ndarray, optional) – square of toroidal Mach-number (Mach-#^2) profile
• m2p (np.ndarray, optional) – derivative of square of Mach-number profile with respect to poloidal flux
• bt2 (np.ndarray, optional) – square of toroidal magnetic field profile
• bp2 (np.ndarray, optional) – square of poloidal magnetic field profile
• b2 (np.ndarray, optional) – total magnetic field (Bt2 + Bp2) profile
• fedge (np.ndarray, optional) – force profile at plasma edge

psi

np.ndarray – poloidal flux

psin

np.ndarray – normalized poloidal flux

ffp

np.ndarray – FF’ profile

f

np.ndarray – F profile

pp

np.ndarray – p and p’ profile

p

np.ndarray – p profile

jb

np.ndarray – <J.B> profile

q

np.ndarray – q profile

rhovn

np.ndarray – normalized rho

epot

np.ndarray – electrostatic potential

m2

np.ndarray – square of toroidal Mach-number (Mach-#^2) profile

m2p

np.ndarray – derivative of square of Mach-number profile with respect to poloidal flux

bt2

np.ndarray – square of toroidal magnetic field profile

bp2

np.ndarray – square of poloidal magnetic field profile

b2

np.ndarray – total magnetic field (Bt2 + Bp2) profile

fedge

np.ndarray – force profile at plasma edge

class lib.lrdfit.StrikePointCoords

Coordinates, derivative, angles, R*Bp and flux expansion of strike points

Strike points are points at where the particles escaping from the plasma and following the field lines outside the seperatrix, collide with the walls.

rvsin

list of float – R coordinates of inner strike points

zvsin

list of float – Z coordinates of inner strike points

rvsout

list of float – R coordinates of outer strike points

zvsout

list of float – Z coordinates of outer strike points

dpsidrin

float – derivative of poloidal flux with respect to R coordinate of inner strike point

dpsidzin

float – derivative of poloidal flux with respect to Z coordinate of inner strike point

dpsidrout

float – derivative of poloidal flux with respect to R coordinate of outer strike point

dpsidzout

float – derivative of poloidal flux with respect to Z coordinate of outer strike point

angpin

float – Angle between plate and Bpol - inner [Degrees]

angpout

float – Angle between plate and Bpol - outer [Degrees]

angtin

float – Angle between plate and Btot - inner [Degrees]

angtout

float – Angle between plate and Btot - outer [Degrees]

rbprmxb

float – R*Bp - Rmax on boundary

rbpin

float – R*Bp - inner

rbpout

float – R*Bp - outer

fexpin

float – Flux expansion - inner

fexpout

float – Flux expansion - outer

compute(x, rb, zb, lim, iplim=None, fedge=None, debug=None)

compute the strike points coordinates

Parameters:

• x (PoloidalFluxProfile) – x-point porfile data
• rb (np.ndarray) – boundary coordinates in R axis
• zb (np.ndarray) – boundary coordinates in Z axis
• lim (PoloidalFluxProfile) – poloidal flux profile of the limiter
• iplim (efit.PoloidalFlux, optional) – plasma current at limiter coordinates
• fedge (np.ndarray, optional) – force at plasma edge (?)
• debug (bool, optional) – option to turn debug mode on or off

class lib.lrdfit.StrikePointCtrl(nspc=0, rsp=None, zsp=None, wpsisp=0.0, wbpsp=0.0)

Strike-point control structure

Parameters:

• nspc (int, optional) – 0 to n to add n boundary-point positions to the shape control algorithm
• rsp (np.ndarray, optional) – R values of boundary-point position
• zsp (np.ndarray, optional) – Z values of boundary-point position
• wpsisp (np.ndarray, optional) – weighting of boundary-point flux
• wbpsp (np.ndarray, optional) – weighting of boundary-point flux gradient

nspc

int – 0 to n to add n boundary-point positions to the shape control algorithm

rsp

np.ndarray – R values of boundary-point position

zsp

np.ndarray – Z values of boundary-point position

wpsisp

np.ndarray – weighting of boundary-point flux

wbpsp

np.ndarray – weighting of boundary-point flux gradient

irsp

np.ndarray – 0 = extrapolate R to nearby limiter position, 1 = use exact R value

izsp

np.ndarray – 0 = extrapolate Z to nearby limiter position, 1 = use exact Z value

imodesp

np.ndarray – 0 = dpsin only, 1 = inverse flux expansion, 2 = B-angle [degrees]

invfexpsp

np.ndarray – value for inverse flux expansion at boundary-point

winvfexpsp

np.ndarray – weighting for inverse flux expansion at boundary-point

bangdegsp

np.ndarray – value of B-field angle of incidence at boundary-point [degrees]

wbangdegsp

np.ndarray – weighting for B-field angle of incidence at boundary-point [degrees]

class lib.lrdfit.TRANSPEqm(trdatread=None, trexpt=None, trshot=None, trtime=None, trdtav=None, trid=None, tryear=None)

TRANSP equilibrium parameter structure

Parameters:

• trdatread (int, optional) – Read TRANSP data?
• trexpt (str, optional) – Read TRANSP data?
• trshot (int, optional) – TRANSP shot number
• trtime (float, optional) – TRANSP shot time duration
• trdtav (float, optional) – TRANSP dtav(?)
• trid (str, optional) – TRANSP ID
• tryear (int, optional) – TRANSP year

trdatread

int – Read TRANSP data?

trexpt

str – Read TRANSP data?

trshot

int – TRANSP shot number

trtime

float – TRANSP shot time duration

trdtav

float – TRANSP dtav(?)

trid

str – TRANSP ID

tryear

int – TRANSP year

Examples

>>> trpars = TRANSPEqm(trdatread=0,
...                    trexpt='NSTX',
...                    trshot=200007,
...                    trtime=0.4,
...                    trdtav=0.05,
...                    trid='Y01',
...                    tryear=2006)

EFIT

class lib.efit.BoundaryShapeSquareness

Geometric moments of the boundary shape and squareness

r0

float – major radii of midplane crossing points

z0

float – vertical center as Z position of Rmax

a

float – Plasma minor radius

ku

float – upper elongation

kl

float – lower elongation

du

float – upper triangularity

dl

float – lower triangularity

squo

np.ndarray – upper outer squareness

squi

np.ndarray – upper inner squareness

sqli

np.ndarray – lower inner squareness

sqlo

np.ndarray – lower outer squareness

xkds

np.ndarray – X coordinate of fitted boundary shape onto reference angle

zkds

np.ndarray – Z coordinate of fitted boundary shape onto reference angle

tkds

np.ndarray – theta of equilibirum boundary

isquse

np.ndarray – indices for computing squareness values

compute(xbnd, zbnd, tfit=None, sqmin=-0.35, sqmax=0.6)

Compute geometric moments of the boundary shape and squareness

Given X and Z on the boundary (xb, zb), compute geometric moments of the boundary shape including squareness using boundary shape (and squareness) definition from: A.D. Turnbull et al., Phys. Plasmas 6, 1113 (1999)

Parameters:

• xbnd (np.ndarray) – X coordinate of boundary
• zbnd (np.ndarray) – Z coordinate of boundary
• tfit (np.ndarray, optional) – theta of equilibirum boundary for fitting
• sqmin (float, optional) – minimum squareness
• sqmax (float, optional) – maximum squareness

class lib.efit.CoordData(shape=None)

Coordinates data in EFIT

Parameters:shape (tuple or int, optional) – shape of coordinate data

x

np.ndarray – data for X coordinate

y

np.ndarray – data for Y coordinate

z

np.ndarray – data for Z coordinate

r

np.ndarray – data for R coordinate

class lib.efit.EqmFluxCoord(nr=101, nt=101, jacobian=None, theta=None, psivec=None, x=None)

Equilibrium solution in flux coordinate

This class contains equilibrium solution in flux coordinate derived from EFIT G-EQDSK data

Parameters:

• nr (int, optional) – number of radius coordinate on surface
• nt (int, optional) – number of theta coordinate on surface
• jacobian (np.ndarray, optional) – jacobian of flux coordinate
• theta (np.ndarray, optional) – uniform poloidal angles in 2pi period
• psivec (np.ndarray, optional) – uniform poloidal flux on grid

coords

str – name of user-defined coordinate system

points

CoordData – points in X, Y, Z, R coordinate

pesttheta

np.ndarray – pest-theta coordinates in the new coordinate system

nt

int – number of poloidal angles

nr

int – number of flux surfaces for toroidal angles

jacs

FluxCoordJacobian – jacobian of flux coordinate

theta

np.ndarray – uniform poloidal angles in 2pi period

psivec

np.ndarray – uniform poloidal flux on grid

time

np.ndarray – time duration array in second

shot

int – shot number

eindex

np.ndarray – EFIT index

boundary

CoordData – boundary values in X, Z coordinate

torflux

np.ndarray – toroidal flux

psi

np.ndarray – poloidal flux

rhopol

np.ndarray – square root of normalized poloidal flux

rhotor

np.ndarray – square root of normalized toroidal flux

psin

np.ndarray – normalized poloidal flux for each poloidal angle

psinh

np.ndarray – “helically distorted” normalized poloidal flux

psinorm

np.ndarray – normalized poloidal flux in flux coordinate

phirad

float – poloidal angle in rad

phideg

float – poloidal angle in degree

Rmin

np.ndarray – minimum points of R coordinate

Rmax

np.ndarray – maximum points of R coordinate

Bmin

np.ndarray – minimum field

Bmax

np.ndarray

epsg

np.ndarray – global inverse aspect ratio profiles

epsl

np.ndarray – local inverse aspect ratio profiles

epsb

np.ndarray – mod-B equivalent inverse aspect ratio profiles

parclen

np.ndarray – poloidal arc-length

derchk

np.ndarray – test of the derivative calculations, which should be unity at (x, z)

p

np.ndarray – p on uniform flux grid

pp

np.ndarray – p’ on uniform flux grid

ffp

np.ndarray – FF’ on uniform flux grid

fp

np.ndarray – F’ on uniform flux grid

f

np.ndarray – F on uniform flux grid

q

np.ndarray – recalculated q profile

q_orig

np.ndarray – original q profile

B

np.ndarray – total magnetic field

Bphi

np.ndarray – toroidal magnetic field

Bpol

np.ndarray – poloidal magnetic field

Br

np.ndarray – R component of magnetic field

Bz

np.ndarray – Z component of magnetic field

saB2

np.ndarray – flux surface average of B^2

saBphi2

np.ndarray – flux surface average of Bphi^2

saBpol2

np.ndarray – flux surface average of Bpol^2

volume

np.ndarray – total volume

rhovol

np.ndarray – square root of normalized volume

dVdpsi

np.ndarray – derivative of volume with respect to poloidal flux

Jchease

np.ndarray – <Jphi/R>/<1/R> for CHEASE option NSTTP=2

Jjsolv

np.ndarray – <J.B>/<B.Grad-phi>(psi) for J-SOLVER

Jphi

np.ndarray – toroidal current density

Jpol

np.ndarray – poloidal current density

saJBfp

np.ndarray – flux surface average of <J.B> (F’ componet)

saJBpp

np.ndarray – flux surface average of <J.B> (p’ componet)

saJBtot

np.ndarray – flux surface average of <J.B> in total

saRm2

np.ndarray – flux surface average of 1/R^2

m2

np.ndarray – toroidal Mach number squared

m2p

np.ndarray – d(M(psi))^2/dpsi

dBphidt

np.ndarray – derivative of toroidal magnetic field with respect to time

dBzdt

np.ndarray – derivative of Z component of magnetic field with respect to time

dBrdt

np.ndarray – derivative of R component of magnetic field with respect to time

vloop

np.ndarray – total loop voltage

vloopt

np.ndarray – toroidal component of loop voltage

vloopp

np.ndarray – poloidal component of loop voltage

sigma

np.ndarray – conductivity

Er

np.ndarray – R component of electric field

Ez

np.ndarray – Z component of electric field

Ephi

np.ndarray – toroidal electric field

Epol

np.ndarray – poloidal electric field

saEBphi

np.ndarray – flux surface averages of E.B (toroidal component)

saEBpol

np.ndarray – flux surface averages of E.B (poloidal component)

saEBtot

np.ndarray – flux surface averages of E.B in total

area_integrate(a, each_flux_surface=False)

Compute area integral

Compute cross-sectional area integral A. dA on each flux surface or from axis to boundary

Parameters:

• a (np.ndarray) – input array A
• each_flux_surface (bool, optional) – compute on each flux surface

Returns:

area integral

Return type:

float or np.ndarray

compute(gs, rexp0=2.0, rexp1=0.25, psifrac=0.997, ctol=0.0001, maxit=10, coords='lr', field_subset=None, linear=False, fast=False, fc=None, qedge=None, pedge=None, rhomin=None, rhomax=None, quiet=None, full_bnd=None, eindex=None)

Compute equilibrium solution

Given an EFIT g-structure, compute the full equilibrium solution in flux coordinates.

Parameters:

• gs (GStructure) – G-EQDSK data
• rexp0 (float, optional) – psi scales as rminor^radexp1 near axis
• rexp1 (float, optional) – psi scales as rminor^radexp0 near edge
• psifrac (float, optional) – fraction of psi spanned by flux surfaces
• ctol (float, optional) – flux surface convergence parameter
• maxit (int, optional) – maximum number of iterations
• coords (str, optional) – default coordinate system
• field_subset (Field, optional) – Magnetic and electric field data
• linear (bool, optional) – if radial coordinate is linear in poloidal flux
• fast (bool, optional) – if show time to calculate magnetic field
• fc (FluxCoord, optional) – flux coordinate data
• qedge (np.ndarray, optional) – specified q-edge to change flux fraction
• pedge (np.ndarray, optional) – force the edge pressure to this value
• rhomin (np.ndarray, optional) – minimum rho
• rhomax (np.ndarray, optional) – maximum rho
• quiet (bool, optional) – if print extra information
• full_bnd (bool, optional) – if use full boundary
• eindex (np.ndarray, optional) – EFIT index

fs_average_array(a)

Compute flux surface average

Compute flux surface average of “A” using Jacobian “J” & angle “theta” on a single flux surface (i.e. A, J, and theta are 1D theta arrays)

Parameters:a (np.ndarray) – input array Returns:flux surface average Return type:np.ndarray

saip_integrate(a, divfree=None)

Compute toroidal plasma current

Given a 2D psi x time array (A) of the surface average <J.B> and the equilibrium data returned by “efitg_flux_coords” (efc), compute the toroidal plasma current enclosed by each flux surface. If keyword “divfree” is set, include only the divergence-free current.

Parameters:

• a (np.ndarray) – input 2D array A
• divfree (bool, optional) – if include only the divergence-free current

Returns:

saIp

Return type:

float

surface_integrate(a, each_flux_surface=False)

Compute surface area integral

Compute surface area integral A.dS on each flux surface or on boundary

Parameters:

• a (np.ndarray) – input array for surface integration
• each_flux_surface (bool, optional) – compute on each flux surface

Returns:

surface integral

Return type:

float

volume_integrate(psi_array, each_flux_surface=False, profile=False, average=False)

Compute volume integral

Given a 2D theta x psi array (A) and the equilibrium data in “EqmFluxCoord”, compute volume integral A*dV

Parameters:

• psi_array (np.ndarray) – 1D or 2D theta x psi array
• each_flux_surface (bool, optional) – compute volume integral on each flux surface
• profile (bool, optional) – Compute volume average of profile f(psin)*dV from axis to boundary
• average (bool, optional) – normalized the volume integral

Returns:

volume integral

Return type:

float or np.ndarray

class lib.efit.Equilibrium

Equilibrium data

converted

bool – if the equilibrium data has been converted from equilibrium solution in flux coordinate.

derived

bool – if the boundary data has been derived from equilibrium data.

x

np.ndarray – X coordinate

z

np.ndarray – Z coordinate

theta

np.ndarray – uniform poloidal angles in 2pi period

psi

np.ndarray – uniform poloidal flux on grid

p

np.ndarray – p on uniform flux grid

pp

np.ndarray – p’ on uniform flux grid

ffp

np.ndarray – FF’ on uniform flux grid

f

np.ndarray – f on uniform flux grid

q

np.ndarray – q on uniform flux grid

m2

np.ndarray – toroidal Mach number squared

m2p

np.ndarray – d(M(psi))^2/dpsi

xbnd

np.ndarray – X coordinate at boundary

zbnd

np.ndarray – Z coordinate at boundary

rhobnd

np.ndarray – rho at boundary

thetabnd

np.ndarray – poloidal angle at boundary

psin

np.ndarray – normalized poloidal flux

rhovol

np.ndarray – square root of normalized volume

rhotor

np.ndarray – square root of normalized toroidal flux

jacobian

np.ndarray – jacobian of flux coordinate

saJdBn

np.ndarray – <J.B> / <B.grad-phi>

saBp2

np.ndarray – flux surface average Bp^2

saRm2

np.ndarray – flux surface average 1/R^2

iJdt

np.ndarray – momentum of jacobian

xaxis

np.ndarray – X coordinate of axis

zaxis

np.ndarray – Z coordinate of axis

r0

float – plasma geometric center

aminor

float – plasma minor radius

volume

float – plasma volume

cxarea

float – plasma cx area

surfarea

float – plasma surface area

lpol

float – boundary poloidal arc length

kappa

float – plasma elongation

kappa0

float – plasma elongation for r –> 0

tritop

float – upper triangularity

tribot

float – lower triangularity

shape_fit

BoundaryShapeSquareness – boundary shape and squareness

ip_pp

float – toroidal current from p’

ip_jb

float – toroidal current from <J.B>

ip

float – total toroidal current

li_ga

float – internal inductance (GA)

bt0

np.ndarray – vacuum toroidal field at R0

betatv0

float – vacuum toroidal beta

betan

float – normalized beta

ctroyon

float – troyon normalized beta

wmhd

float – 3/2 * volume integral of plasma pressure

ptap

float – peak to average pressure

q0

float – q at magnetic axis

qmin

float – minimum q(psi)

rhoqmin

float – square root of psin of minimum q

q95

float – safety factor at psin = 0.95

q99

float – safety factor at psin = 0.99

qstar

float – q*

rhoq

np.ndarray – SQRT(psin) of integer q values 1.0 - 5.0

dqdr_q

np.ndarray – dq/SQRT(psin) atinteger q values 1.0 - 5.0; dqdr_q1 - dpdr_q5 in original IDL.

dqdr_r80

float – derivative of q with respect to rho for rhoh <= 0.8; rhoh is square root of “helically distorted” normalized flux.

dqdr_r90

float – derivative of q with respect to rho for rho <= 0.9

q_r80

float – q for rhoh = 0.8

q_r90

float – q for rhoh = 0.9

psitor

np.ndarray – enclosed actual toroidal flux

psitorv

np.ndarray – enclosed vacuum toroidal flux

fdiamag

np.ndarray – diamagnetic flux

convert_efc_to_eqs(efc)

Convert equilibrium solution from flux coordinate to compute boundary and other data

Parameters:efc (EqmFluxCoord) – equilibrium solution in flux coordinate converting to equilibrium data

derive(shape_fit=False)

Derive boundary and other data from equilibrium data

Refer to “eqs_derived_data()” from IDL library “efit_routines_jem.pro”

Parameters:shape_fit (BoundaryShapeSquareness, optional) – Geometric moments of the boundary shape and squareness

class lib.efit.Field

Electric and magnetic field data

psi

np.ndarray – poloidal flux

Bphi

np.ndarray – toroidal magnetic field

dpsidt

np.ndarray – derivative of poloidal flux with respect to time

dBphidt

np.ndarray – derivative of toroidal magnetic field with respect to time

Er

np.ndarray – R component of electric field

Ez

np.ndarray – Z component of electric field

points

CoordData – R, Z coordinate

class lib.efit.FluxCoord(nradius, ntheta)

Flux coordinate data

Flux coordinates [2] in the context of magnetic confinement fusion (MCF) is a set of coordinate functions adapted to the shape of the flux surfaces of the confining magnetic trap. They consist of one flux label, often termed psi and two angle-like variables theta, psi whose constant contours on the flux (psi(x)=constant) surfaces close either poloidaly (theta) or toroidallly (phi).

In this coordinates, equilibrium vector fields like the magnetic field B or current density j have simplified expressions. A particular kind of flux coordinates, generally called magnetic coordinates, simplify the B-field expression further by making field lines look straight in the (theta, phi) plane of that family of coordinates.

Parameters:

• nradius (int) – number of radius points on surface
• ntheta (int) – number of theta points on surface

nrad

int – number of radius points on surface

nthe

int – number of theta points on surface

xfit

np.ndarray – X coordinate of the fitting grid

zfit

np.ndarray – Z coordinate of the fitting grid

psifit

np.ndarray – fitted poloidal flux on the fitting grid

psinorm

np.ndarray – normalized poloidal flux on the fitting grid

xbnd

np.ndarray – X coordinate of boundary for the new surface

zbnd

np.ndarray – Z coordinate of boundary for the new surface

angle

np.ndarray – poloidal angle of x,z coordinates

theta

np.ndarray – interpolated theta for [0, 2*pi] interval

qoverf

np.ndarray – q / F where F = R * Bt

References

[2]	Fusion Wiki!.

change_theta_coordinate(thev, ntheta=None, coords=None)

Change the coordinate of theta

Given the 2D arrays X & Z and 1D arrays theta and psi, re-map the X,Z arrays in the theta variable - possibly with a different Jacobian. METHOD: Conservation of volume ==> J1 * dthe1 = J2 * dthe2(the1). Since J2 is specified, dthe2 = J1 / J2 * dthe1. Then, compute anti-derivative of dthe2 and re-normalize.

Parameters:

• thev (np.ndarray) – theta variable to remap X, Z arrays
• ntheta (int, optional) – number of points in the new theta coordinate.
• coords (str, optional) –

  new coordinate system:

  ’pt’ –> “pest theta” poloidal coordinates

  ’ea’ –> “equal-arc” poloidal coordinates

  ’ha’ –> “hamada” poloidal coordinates

  ’lr’ –> “linear-ray” poloidal coordinates

compute(psidata, r_bnd, z_bnd, r_axis, z_axis, radexp0, radexp1, fluxfrac, rscale, dconv, miters, coords='lr', rhomin=None, rhomax=None, full_bnd=False)

Compute flux coordinate based on poloidal flux data

Given the reformed poloidal flux data, boundary and axis info, and interpolation parameters, return flux coordinates X & Z (and a bunch of other stuff)

Parameters:

• psidata (PoloidalFlux) – poloidal flux data
• r_bnd (np.ndarray) – vector of boundary major radius points
• z_bnd (np.ndarray) – vector of boundary Z points
• r_axis (np.ndarray) – major radius of magnetic axis
• z_axis (np.ndarray) – vertical position of magnetic axis
• radexp0 (float) – psi scales as rminor^radexp0 near axis
• radexp1 (float) – psi scales as rminor^radexp1 near edge
• fluxfrac (float) – fraction of total flux spanned by surfaces
• rscale (float) – extrapolate past boundary this far when interpolating
• dconv (float) – flux surface convergence parameter
• miters (int) – maximum number of iterations
• coords (str, optional) – coordinate system choices
• rhomin (np.ndarray, optional) – minimum rho
• rhomax (np.ndarray, optional) – maximum rho
• full_bnd (bool, optional) – if use full boundary

fit_boundary(xb, zb, udsym=False, uniform=False)

Fit new boundary data

Given raw boundary information from EFIT, re-order the boundary points and interpolate in theta

Parameters:

• xb (np.ndarray) – raw boundary data in X axis
• zb (np.ndarray) – raw boundary data in Z axis
• udsym (bool, optional) – if symmetrically fit?
• uniform (bool, optional) – if uniformly fit?

class lib.efit.FluxCoordJacobian

Jacobian for flux coordinate

jacobian

np.ndarray – jacobian matrix - 2D array

xthe

np.ndarray – derivative of X with respect to poloidal angle

xpsi

np.ndarray – derivative of X with respect to poloidal flux

rthe

np.ndarray – derivative of X with respect to poloidal angle

rpsi

np.ndarray – derivative of X with respect to poloidal flux

zthe

np.ndarray – derivative of Z with respect to poloidal angle

zpsi

np.ndarray – derivative of Z with respect to poloidal flux

delstar_psi

np.ndarray – rho derivative for delstar-psi

moments

JacobianMoments – moments of the Jacobian

compute(x, z, thev, psiv, delstar_psi=None)

Compute jacobian for flux coordinate

Given the 2D arrays X & Z and 1D arrays theta and psi, compute the flux coordinate Jacobian and related derivatives on the X,Z grid. Arrays must repeat themselves in the theta direction spanning [0,2PI] (i.e. strictly periodic). Del-star psi on the grid is optionally computed with the “delstar_psi” keyword.

Parameters:

• x (np.ndarray) – 2D array X coordinate
• z (np.ndarray) – 2D array Z coordinate
• thev (np.ndarray) – 1D array data theta variable
• psiv (np.ndarray) – 1D array data psi variable
• delstar_psi (bool, optional) – if compute del-star psi data?

class lib.efit.GFileName

G-EQDSK file name

shotstr

str – string contains shot number

timestr

str – string contains time

timestr_

str – string contains full length time

gname

str – file name for G-EQDSK format data

gname_

str – G-EQDSK file name with full length time string

generate(shot, time, afile=False)

Generate file name

Parameters:

• shot (np.ndarray) – shot number
• time (np.ndarray) – time value
• afile (bool, optional) – if use ‘a’ as prefix

class lib.efit.GStructure

G-EQDSK structure for poloidal flux data

A standard format for storing tokamak equilibrium data is G-EQDSK [3] which contains the poloidal flux in (R,Z) and 1D profiles of pressure, f = R * B(phi), safety factor q, and other quantities related to the solution.

source

str – path of source data file

shot

int – shot number

time

np.ndarray – time base to use for bound determination

ecase

list of np.ndarray – identification character string

error

np.ndarray – error value

mw

np.ndarray – number of horizontal R grid points

mh

np.ndarray – number of vertical Z grid points

xdim

np.ndarray – horizontal dimension in meter of computational box

zdim

np.ndarray – vertical dimension in meter of computational box

rzero

np.ndarray – initial R value

rgrid1

np.ndarray – R of left side of domain

zmid

np.ndarray – Z of center of computational box in meter

rmaxis

np.ndarray – R of magnetic axis in meter

zmaxis

np.ndarray – Z of magnetic axis in meter

ssimag

np.ndarray – poloidal flux at magnetic axis in Weber /rad

ssibry

np.ndarray – poloidal flux at the plasma boundary in Weber /rad

bcentr

np.ndarray – vacuum toroidal magnetic field in Tesla at RCENTR

cpasma

np.ndarray – plasma current

fpol

np.ndarray – poloidal current function in m-T, F = RBT on flux grid

pres

np.ndarray – Plasma pressure in nt / m^2 on uniform flux grid

ffprim

np.ndarray – FF’(psi) in (mT)^2 (Weber/rad) on uniform flux grid

pprime

np.ndarray – P’(psi) in (nt/m^2) / (Weber/rad) on uniform flux grid

psirz

np.ndarray – poloidal flux in Weber / rad on the rectangular grid points

qpsi

np.ndarray – q values on uniform flux grid from axis to boundary

nbdry

np.ndarray – number of boundary points

limitr

np.ndarray – number of limiter points

bdry

np.ndarray – R, Z coordinate of boundary points

lim

np.ndarray – X, Y coordinate of limiter

r

np.ndarray – radial coordinates of grid

z

np.ndarray – elevation coordinates of grid

rhovn

np.ndarray – normalized rho

epoten

np.ndarray – electrostatic potential function Er/RBz

m2

np.ndarray – square of toroidal Mach-number (Mach-#^2) profile

m2p

np.ndarray – derivative of square of Mach-number profile with respect to poloidal flux

pvsrz0

np.ndarray – variables for dynamic pressure model

pwvsrz0

np.ndarray – variables for dynamic pressure model

ptvsrz0

np.ndarray – variables for dynamic pressure model

pdvsrz0

np.ndarray – variables for dynamic pressure model

qvsrz0

np.ndarray – variables for dynamic pressure model

ptrz

np.ndarray – variables for dynamic pressure model

References

[3]	G-EQDSK file format documentation.

mds_gstructures(shot, tree='efit', index=1, tmin=None, tmax=None, twant=None, plasma=False, auto=False, nw=None, nh=None)

Read G-EQDSK data from MDS+ tree

This function reads data from the EFIT MDS+ tree. However, the dimensions are not hard-wired and are generally incompatible with those required by readg().

NOTE 1: tmin,tmax above have units of seconds

NOTE 2: If tmax <= tmin, time-slice closest to tmin is returned

Parameters:

• shot (int) – shot number
• tree (str, optional) – MDS EFIT tree name, default is ‘efit’
• index (int, optional) – child tree index, default is 1
• tmin (float, optional) – minimum time of returned data default is MIN(time)
• tmax (float, optional) – maximum time of returned data default is MAX(time)
• twant (np.ndarray, optional) – if set, return arrays with times closest to elements of keyword “twant”
• plasma (bool, optional) – if set, return only time-slices with closed flux
• auto (bool, optional) – see auto keyword for rescale_gstructures
• nw (np.ndarray, optional) – use interpolation to redimension mw –> nw
• nh (np.ndarray, optional) – use interpolation to redimension mh –> nh

readg_eqfit(filepath)

Read G-EQDSK file

The method is ported from readg_lrdfit.pro: readg_file This method read in local data file in G-EQDSK format

Parameters:filepath (str) – file path to data

redimension_gstructures(nw, nh)

Change dimension of G-EQDSK data structure

Use interpolation to re-dimension a G-EQDSK structure.

Parameters:

• nw (int) – new dimension in radial and poloidal flux
• nh (int) – new dimension in vertical/height

rescale_gstructures(ip_sign=None, auto=False, qpositive=False, qscale=None, cpasma_positive=False, bt_sign=None, psiswap=False, rpshift=None, zpshift=None, rpscale=None, zpscale=None, limstruc=None, bscale=None)

Rescale G-EQDSK data

This function re-scales the poloidal flux to be consistent with the specified sign of Ip, where Ip is the toroidal component of the plasma current in cylindrical coordinates, i.e. + –> C.C.W. as viewed from above.

Parameters:

• ip_sign (float, optional) – sing of Ip
• auto (bool, optional) – If the “auto” keyword is specified, it is assumed that “cpasma” has the correct sign, the poloidal flux is re-scaled accordingly, and q is forced to be positive.
• qpositive (bool, optional) – Note that in the re-definition of q, the theta flux coordinate is C.C.W. in the poloidal plane, so that if Ip > 0 and Bt > 0, q will be negative, since then B.grad(theta) < 0. This effect is bypassed if the “qpositive” keyword is used which forces q to be positive.
• qscale (float, optional) – If keyword “qscale” is defined, the q profile is scaled by the factor qscale while still satisfying the equilibrium constraint. Typically, qscale controls q(edge) for paramagnetic plasmas and q(0) for diagmagnetic plasmas
• cpasma_positive (bool, optional) – if cpasma is positive?
• bt_sign (float, optional) – sign of toroidal field
• psiswap (bool, optional) – Swap sign of poloidal flux and corresponding derivatives
• rpshift (np.ndarray, optional) – shift plasma R by this amount
• zpshift (np.ndarray, optional) – shift plasma Z by this amount
• rpscale (np.ndarray, optional) – scale plasma R by this amount
• zpscale (np.ndarray, optional) – scale plasma Z by this amount
• limstruc (optional) – use different limiter contains limstruc.rlim and limstruc.zlim
• bscale (np.ndarray, optional) – scale Bt (and F, Ip, psi, p) conserving q and beta

write_gfile(path=None, extension=None, eqname=None, gz=False, rotation=False, microsec=False, namefile=None)

Write G-EQDSK file

This procedure takes a single g-structure and generates an ASCII geqdsk file which can be read by the readg() function above. This routine provides backward compatibility for codes which use ASCII g-files as input and can’t read the MDS tree directly.

Parameters:

• path (str, optional) – file path
• extension (str, optional) – file extension
• eqname (str, optional) – time string
• gz (bool, optional) – if the file is packed as gz-file
• rotation (bool, optional) – if write rotation data
• microsec (bool, optional) – if use micro second as time unit
• namefile (str, optional) – user-defined file name

class lib.efit.JacobianMoments(npsi=None)

Moments of the Jacobian

Parameters:npsi (int, optional) – number of coordinates for poloidal flux

iJdt

np.ndarray

iJBp2dt

np.ndarray

iJRm2dt

np.ndarray

class lib.efit.PoloidalFlux(psi=None, r=None, z=None, nr=None, nz=None)

Poloidal flux data

Parameters:

• psi (np.ndarray, optional) – poloidal flux
• r (np.ndarray, optional) – R coordinate of poloidal flux
• z (np.ndarray, optional) – Z coordinate of poloidal flux
• nr (int, optional) – number of R coordinate
• nz (int, optional) – number of Z coordinate

psi

np.ndarray – poloidal flux

r

np.ndarray – R coordinate of poloidal flux

z

np.ndarray – Z coordinate of poloidal flux

nr

int – number of R coordinate

nz

int – number of Z coordinate

dr

np.ndarray – derivative of R coordinate

dz

np.ndarray – derivative of Z coordinate

rmin

np.ndarray – minimum R coordinate

zmin

np.ndarray – minimum Z coordinate

dpsidr

np.ndarray – derivative of poloidal flux with respect to R

dpsidz

np.ndarray – derivative of poloidal flux with respect to Z

d2psidrdz

np.ndarray – second derivative of poloidal flux with respect to R and Z

ixpoint

bool – if found possible poloidal field X point at R, Z (m)

iopoint

bool – if found possible poloidal field O point at R, Z (m)

gradpsi()

Compute gradient of poloidal flux

computes grad-psi on solution grid for subsequent 2D interpolation, if needed This method is not optmized for performance

interpolate_psi(newr, newz)

Compute the interpolated flux function at new R, Z

Parameters:

• newr (np.ndarray) – new R for interpolation
• newz (np.ndarray) – new Z for interpolation

Returns:

a new instance of class PoloidalFlux with interpolated data

Return type:

PoloidalFlux

process(gs)

Convert g-structure data into poloidal flux data

The following routines are for general post-processing of g-structure data. They are used primarily for computing which array indices are enclosed by the plasma bounary and for refining the boundary data. Reforms the raw flux and R,Z grids from EFIT

Parameters:gs (GStructure) – raw g-structure data

refine_null_position(rin, zin, quiet=False)

Refine null point position

Refines the position of a poloidal field null point using bi-cubic interpolation Two methods are used simulaneously to find both X-points and O-points

Parameters:

• rin (np.ndarray) – input r axis coordinate
• zin (np.ndarray) – input z axis coordinate
• quiet (bool, optional) – option to control command line output

Returns:

• np.ndarray – R coordinate of found null points
• np.ndarray – Z coordinate of found null points

class lib.efit.RotationModifiedPressure

Pressure data modified by rotation

p_psi_R

np.ndarray – 2D pressure including density variation = p(psi,R) on [theta,rho] grid

pp_psi_R

np.ndarray – 2D dp(psi,R)/dpsi at fixed R

sa_p_psi_R

np.ndarray – Flux surface averages of p(psi,R) at fixed R

sa_pp_psi_R

np.ndarray – Flux surface averages of dp(psi,R)/dpsi at fixed R

rotate(p_psi, pp_psi, m2_psi, m2p_psi, r_2d, internal_efc)

Rotate and modify the pressure data

Parameters:

• p_psi (np.ndarray) – pressure as a function of poloidal flux
• pp_psi (np.ndarray) – pressure prime as a function of poloidal flux
• m2_psi (np.ndarray) – square of toroidal Mach-number as a function of poloidal flux
• m2p_psi (np.ndarray) – derivative of square of Mach-number profile with respect to poloidal flux
• r_2d (np.ndarray) – 2D R array
• internal_efc (EqmFluxCoord) – internal subset of equilibrium solution in flux coordinates

class lib.efit.SurfaceComparisonParams(refsurf=None)

Parameters for surface comparison

Parameters:refsurf (CoordData) – coordinate of reference surface shape

rmsdev

np.ndarray – square root of mean of surface matching errors

absdev

np.ndarray – mean of square root of surface matching errors

xref

np.ndarray – X coordinate of reference surface shape

zref

np.ndarray – Z coordinate of reference surface shape

xcom

np.ndarray – X coordinate of comparison surface shape

zcom

np.ndarray – Z coordinate of comparison surface shape

compute(comsurf)

Compute deviation of reference surface shape

Given a reference surface shape (refsurf), compute the deviation of that shape from a “comparison” surface (comsurf) and re-map the comparison surface to have the same poloidal angle distribution as the reference surface. The surfaces can have different numbers of elements and poloidal coordinates.

Parameters:comsurf (CoordData) – coordinate of surface shape for comparison

NEOC

class lib.neoc.BootstrapCurrent

Collisionless bootstrap current

This library is mostly derived from following article

[4] Neoclassical conductivity and bootstrap current formulas for general axisymmetric equilibria and arbitrary collisionality regime

Physics of Plasmas 6, 2834 (1999)

15. Sauter and C. Angioni and Y. R. Lin-Liu

Jbscoll

np.ndarray – <J.B> / <B.grad-phi>

Jbscles

np.ndarray – collisionless form of bootstrap current for comparison

Jinduct

np.ndarray – inductive current density

Vloop

np.ndarray – loop voltage

Jeq

np.ndarray – <J.B>/<B.Grad-phi>(psi) for J-SOLVER

ip_total

np.ndarray – total plasma current

ip_eqj

np.ndarray – plasma current from equivalent J

ip_gradp

np.ndarray – Diamag + PS current

ip_bscoll

np.ndarray – Bootstrap current (coll)

ip_bscles

np.ndarray – Bootstrap current (cles)

ip_induct

np.ndarray – inductive current

thprofs

ThermalProfile – thermal profile data

References

[4]	Neoclassical conductivity and bootstrap….

nc_currents_efc(efc, te, pe, ti, pi, zave, zeff, vloop, nc=None, multi=None)

Derive neoclassical bootstrap current from equilibrium solution in flux coordinate

Parameters:

• efc (EqmFluxCoord) – equilibrium solution in flux coordinate
• te (np.ndarray) – electron temperature
• pe (np.ndarray) – electron pressure
• ti (np.ndarray) – ion temperature
• pi (np.ndarray) – ion pressure
• zave (np.ndarray) – average ion charge
• zeff (np.ndarray) – effective ion charge
• vloop (np.ndarray) – loop voltage
• nc (bool, optional) – if derive from neoclassical data
• multi (PlasmaParams, optional) – plasma parameters for multi-ions

class lib.neoc.EquilibriumEpsilon

Epsilon of equilibrium

aminor

np.ndarray – plasma minor radius

Rmajor

np.ndarray – R coordinate of plasma major radius

epsilon

np.ndarray – ratio between minor and major radius

q

np.ndarray – q profile at boundary

compute(fc)

Compute epsilon of equilibrium

Parameters:fc (EqmFluxCoord) – equilibrium solution in flux coordinate

class lib.neoc.NeocConductivity

Neoclassical conductivity

sigma_nc_h

np.ndarray – neoclassical conductivity for h

sigma_nc_s

np.ndarray – neoclassical conductivity for s

sigma_spitzer

np.ndarray – neoclassical Spitzer conductivity

nuestar_s

np.ndarray – arbitrary collisionality nu_e_* for s

nuistar_s

np.ndarray – arbitrary collisionality nu_i_* for s

nuestar_h

np.ndarray – arbitrary collisionality nu_e_* for h

ft

np.ndarray – trapped particle fraction

fteff_h

np.ndarray – effective trapped particle fraction for h

fteff_s

np.ndarray – effective trapped particle fraction for s

epsilon

np.ndarray – ratio between minor and major radius

compute(fc, edens, te, z)

Parameters:

• fc (EqmFluxCoord) – equilibrium solution in flux coordinate
• edens (float or np.ndarray) – electron density [m^-3]
• te (float or np.ndarray) – electron temperature [eV]
• z (float or np.ndarray) – ion charge

class lib.neoc.PlasmaParams(nr=None, ns=None, n=None, t=None, z=None, m=None)

Plasma parameters

nr

int

ns

int

n

float or np.ndarray – density (m^-3)

T

float or np.ndarray – temperature (eV)

Z

float or np.ndarray – ion charge

m

float or np.ndarray – ion mass relative to proton mass

log_lambda_e()

Derive log(lambda) from electron parameters

Returns:log(lambda) derived from electron parameters Return type:float or np.ndarray

log_lambda_i()

Derive log(lambda) from ion parameters

Returns:log(lambda) derived from electron parameters Return type:float or np.ndarray

log_lambda_ii(n2x, t2x, z2, m2)

Derive log(lambda) from two ions

Multi-ion version - from NRL handbook [5]

Parameters:

• n2x (float or np.ndarray) – another ion density [m^-3]
• t2x (float or np.ndarray) – another ion temperature [eV]
• z2 (float or np.ndarray) – another ion charge
• m2 (float or np.ndarray) – another ion mass (relative to proton mass)

Returns:

log(lambda) from multiple ions

Return type:

float or np.ndarray

References

[5]	NRL Plasma Formulary.

class lib.neoc.ThermalProfile

Thermal profile data

psin

np.ndarray – normalized poloidal flux

rhopol

np.ndarray – square root of normalized poloidal flux

rhotor

np.ndarray – square root of normalized toroidal flux

pth

np.ndarray – thermal pressure

pe

np.ndarray – electron pressure

pi

np.ndarray – ion pressure

Te

np.ndarray – electron temperature

Ti

np.ndarray – ion temperature

Zave

np.ndarray – ionization of nuclei –> average ionization [7]

Zeff

np.ndarray – effective ion charge

ftrap

np.ndarray – trapped fraction

nue

np.ndarray – arbitrary collisionality nu_e

nui

np.ndarray – arbitrary collisionality nu_i

cloge

np.ndarray – log(lambda) derived from electron parameters

clogi

np.ndarray – log(lambda) derived from ion parameters

eta_snc

np.ndarray – neoclassical Spitzer resistivity for s

eta_hnc

np.ndarray – neoclassical Spitzer resistivity for h

eta_sp

np.ndarray – Spitzer resistivity

References

[7]	Erratum PoP Vol. 9, No 12 (2002) 5140.

class lib.neoc.TrappedFraction

Trapped particle fraction

ft

np.ndarray – trapped particle fraction

ftu

np.ndarray – upper bound of trapped particle fraction

ftl

np.ndarray – lower bound of trapped particle fraction

compute(flux_coords)

Compute trapped particle fraction using method of Lin-Liu and Miller [6]

“Upper and lower bounds of the effective trapped particle fraction in general tokamak equilibria”

Physics of Plasmas 2, 1666 (1998)

25. 18. Lin-Liu and R. L. Miller

Parameters:flux_coords (EqmFluxCoord) – flux coordinate data

References

[6]	Physics of Plasmas 2, 1666 (1995).

EQSC

class lib.eqsc.Equilibrium(eqs=None)

Equilibrium solution data

Parameters:eqs (efit.Equilibrium, optional) – use poloidal flux and coordinate data from exisiting equilibrium solution to initialize current instance

get_jsolver_eqdata(eqname, ipsign=1.0, btsign=1.0)

Get equilibrium solution from J-solver

Parameters:

• eqname (str) – name of J-solver data file
• ipsign (float, optional) – sign of plasma current
• btsign (float, optional) – sign of toroidal field

psi_theta_interp(xwant, zwant)

Interpolate poloidal angle and flux data

originally ported from IDL source code file: eqsc_routines_jem.pro Using the data returned from get_jsolver_eqdata (eqs), update the flux and angle (using simple interpolation) at the desired rectilinear coordinates (x,z) = (xwant,zwant).

Parameters:

• xwant (np.ndarray) – desired rectilinear coordinates x
• zwant (np.ndarray) – desired rectilinear coordinates z

Indices and tables

• Index
• Module Index
• Search Page