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:

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.

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.

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)

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)

1. 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

PyISOLVER

Overview

How to Acquire PyISOLVER

How to Run

How to Change Input Parameters

How to Plot Figures

Reference to Original IDL ISOLVER

Example Run

Performance

Documentation for the Code

LRDFIT

EFIT

NEOC

EQSC

Indices and tables

Table Of Contents

Related Topics

This Page