Aurora modules¶

Submodules¶

aurora.core module¶

This module includes the core class to set up simulations with aurora. The aurora_sim takes as input a namelist dictionary and a g-file dictionary (and possibly other optional argument) and allows creation of grids, interpolation of atomic rates and other steps before running the forward model.

class aurora.core.aurora_sim(namelist, geqdsk=None, nbi_cxr=None)[source]¶

Bases: object

Class to setup and run aurora simulations.

calc_Zeff()[source]¶

Compute Zeff from each charge state density, using the result of an Aurora simulation. The total Zeff change over time and space due to the simulated impurity can be simply obtained by summing over charge states

Results are stored as an attribute of the simulation object instance.

centrifugal_asym(omega, Zeff, plot=False)[source]¶

Estimate impurity poloidal asymmetry effects from centrifugal forces. See notes the centrifugal_asym() function docstring for details.

In this function, we use the average Z of the impurity species in the Aurora simulation result, using only the last time slice to calculate fractional abundances. The CF lambda factor

Parameters

omega – array (nt,nr) or (nr,) [ rad/s ] Toroidal rotation on Aurora temporal time_grid and radial rhop_grid (or, equivalently, rvol_grid) grids.
Zeff – array (nt,nr), (nr,) or float Effective plasma charge on Aurora temporal time_grid and radial rhop_grid (or, equivalently, rvol_grid) grids. Alternatively, users may give Zeff as a float (taken constant over time and space).

Keyword Arguments

plot – bool If True, plot asymmetry factor $\lambda$ vs. radius

Returns

array (nr,): Asymmetry factor, defined as $\lambda$ in the centrifugal_asym() function docstring.

Return type

CF_lambda

check_conservation(plot=True, axs=None, plot_resolutions=False)[source]¶

Check particle conservation for an aurora simulation.

Args :

plotbool, optional: If True, plot time histories in each particle reservoir and display quality of particle conservation.
axsmatplotlib.Axes instances, optional: Axes to pass to check_particle_conserv() These may be the axes returned from a previous call to this function, to overlap results for different runs.

Returns :

outdict: Dictionary containing density of particles in each reservoir.
axsmatplotlib.Axes instances , only returned if plot=True: New or updated axes returned by check_particle_conserv()

get_aurora_kin_profs(min_T=1.01, min_ne=10000000000.0)[source]¶: Get kinetic profiles on radial and time grids.

get_par_loss_rate(trust_SOL_Ti=False)[source]¶

Calculate the parallel loss frequency on the radial and temporal grids [1/s].

trust_SOL_Ti should generally be set to False, unless specific Ti measurements are available in the SOL.

get_time_dept_atomic_rates()[source]¶: Obtain time-dependent ionization and recombination rates for a simulation run. If kinetic profiles are given as time-independent, atomic rates for each time slice will be set to be the same.

interp_kin_prof(prof)[source]¶: Interpolate the given kinetic profile on the radial and temporal grids [units of s]. This function extrapolates in the SOL based on input options using the same methods as in STRAHL.

plot_resolutions()[source]¶: Convenience function to show time and spatial resolution in Aurora simulation setup.

run_aurora(D_z, V_z, times_DV=None, nz_init=None, alg_opt=1, evolneut=False, use_julia=False, plot=False)[source]¶

Run a simulation using inputs in the given dictionary and D,v profiles as a function of space, time and potentially also ionization state. Users may give an initial state of each ion charge state as an input.

Results can be conveniently visualized with time-slider using

aurora.slider_plot(rhop,time, nz.transpose(1,2,0),
                   xlabel=r'$\rho_p$', ylabel='time [s]',
                   zlabel=r'$n_z$ [cm$^{-3}$]', plot_sum=True,
                   labels=[f'Ca$^{{{str(i)}}}$' for i in np.arange(nz_w.shape[1]])

Parameters

D_z – arrays, shape of (space,time,nZ) or (space,time) or (space,) Diffusion and convection coefficients, in units of cm^2/s and cm/s, respectively. This may be given as a function of (space,time) or (space,nZ, time), where nZ indicates the number of charge states. If D_z and V_z are found to be have only 2 dimensions, it is assumed that all charge states should have the same transport coefficients. If they are only 1-D, it is further assumed that they are time-independent. Note that it is assumed that D_z and V_z profiles are already on the self.rvol_grid radial grid.
V_z – arrays, shape of (space,time,nZ) or (space,time) or (space,) Diffusion and convection coefficients, in units of cm^2/s and cm/s, respectively. This may be given as a function of (space,time) or (space,nZ, time), where nZ indicates the number of charge states. If D_z and V_z are found to be have only 2 dimensions, it is assumed that all charge states should have the same transport coefficients. If they are only 1-D, it is further assumed that they are time-independent. Note that it is assumed that D_z and V_z profiles are already on the self.rvol_grid radial grid.

Keyword Arguments

times_DV – 1D array, optional Array of times at which D_z and V_z profiles are given. By Default, this is None, which implies that D_z and V_z are time independent.
nz_init – array, shape of (space, nZ) Impurity charge states at the initial time of the simulation. If left to None, this is internally set to an array of 0’s.
alg_opt – int, optional If alg_opt=1, use the finite-volume algorithm proposed by Linder et al. NF 2020. If alg_opt=1, use the older finite-differences algorithm in the 2018 version of STRAHL.
evolneut – bool, optional If True, evolve neutral impurities based on their D,V coefficients. Default is False, in which case neutrals are only taken as a source and those that are not ionized immediately after injection are neglected. This option is NOT CURRENTLY RECOMMENDED, because this method is still under development/ examination.
use_julia – bool, optional If True, run the Julia pre-compiled version of the code. Run the julia makefile option to set this up. Default is False (still under development)
plot – bool, optional If True, plot density for each charge state using a convenient slides over time and check particle conservation in each particle reservoir.

Returns

array, (nr,nZ,nt): Charge state densities [:math::cm^{-3}] over the space and time grids.
N_wallarray (nt,): Number of particles at the wall reservoir over time.
N_divarray (nt,): Number of particles in the divertor reservoir over time.
N_pumparray (nt,): Number of particles in the pump reservoir over time.
N_retarray (nt,): Number of particles temporarily held in the wall reservoirs.
N_tsuarray (nt,): Edge particle loss [:math::cm^{-3}]
N_dsuarray (nt,): Parallel particle loss [:math::cm^{-3}]
N_dsularray (nt,): Parallel particle loss at the limiter [:math::cm^{-3}]
rcld_ratearray (nt,): Recycling from the divertor [:math::s^{-1} cm^{-3}]
rclw_ratearray (nt,): Recycling from the wall [:math::s^{-1} cm^{-3}]

Return type

nz

setup_grids()[source]¶: Method to set up radial and temporal grids given namelist inputs.

setup_kin_profs_depts()[source]¶: Method to set up Aurora inputs related to the kinetic background from namelist inputs.

aurora.atomic module¶

Collection of classes and functions for loading, interpolation and processing of atomic data. Refer also to the adas_files.py script.

class aurora.atomic.CartesianGrid(grids, values)[source]¶

Bases: object

Linear multivariate Cartesian grid interpolation in arbitrary dimensions This is a regular grid with equal spacing.

class aurora.atomic.adas_file(filepath)[source]¶

Bases: object

Read ADAS file in ADF11 format over the given density and temperature grids. Note that such grids vary between files, and the species they refer to may too.

Refer to ADAS documentation for details on each file.

load()[source]¶

plot(fig=None, axes=None)[source]¶

aurora.atomic.balance(logTe_val, cs, n0_by_ne, logTe_, S, R, cx)[source]¶: Evaluate balance of effective ionization, recombination and charge exchange at a given temperature.

aurora.atomic.get_adas_file_types()[source]¶

Obtain a description of each ADAS file type and its meaning in the context of Aurora. For background, refer to:

Summers et al., "Ionization state, excited populations and emission of impurities
in dynamic finite density plasmas: I. The generalized collisional-radiative model for
light elements", Plasma Physics and Controlled Fusion, 48:2, 2006

Returns: Dictionary with keys given by the ADAS file types and values giving a description for them.

aurora.atomic.get_atom_data(imp, filetypes=['acd', 'scd'], filenames=[])[source]¶

Collect atomic data for a given impurity from all types of ADAS files available or for only those requested.

Parameters

imp – str Atomic symbol of impurity ion.
filetypes – list or array-like ADAS file types to be fetched. Default is [“acd”,”scd”] for effective ionization and recombination rates (excluding CX).
filenames – list or array-like, optional ADAS file names to be used in place of the defaults given by adas_file_dict(). If left empty, such defaults are used. Note that the order of filenames must be the same as the one in the “filetypes” list.

Returns

dict: Dictionary containing data for each of the requested files. Each entry of the dictionary gives log-10 of ne, log-10 of Te and log-10 of the data as attributes atom_data[key].logNe, atom_data[key].logT, atom_data[key].data

Return type

atom_data

aurora.atomic.get_cooling_factors(atom_data, logTe_prof, fz, plot=True, ax=None)[source]¶

Calculate cooling coefficients for the given fractional abundances and kinetic profiles.

Parameters

atom_data – dict Dictionary containing atomic data as output by get_atom_data() for the atomic processes of interest. “prs”,”pls”,”plt” and “prb” are required by this function.
logTe_prof – array (nt,nr) Log-10 of electron temperature profile (in eV)
fz – array (nt,nr) Fractional abundances for all charge states of the ion of “atom_data”
plot – bool If True, plot all radiation components, summed over charge states.
ax – matplotlib.Axes instance If provided, plot results on these axes.

Returns

array (nt,nr): Line radiation in the SXR range for each charge state
prsarray (nt,nr): Continuum radiation in the SXR range for each charge state
plttarray (nt,nr): Line radiation (unfiltered) for each charge state. NB: this corresponds to the ADAS “plt” files. An additional “t” is added to the name to avoid conflict with the common matplotlib.pyplot short form “plt”
prbarray (nt,nr): Continuum radiation (unfiltered) for each charge state

Return type

pls

aurora.atomic.get_cs_balance_terms(atom_data, ne_cm3=50000000000000.0, Te_eV=None, maxTe=10000.0, include_cx=True)[source]¶

Get S, R and cx on the same logTe grid.

Parameters

atom_data – dictionary of atomic ADAS files (only acd, scd are required; ccd is necessary only if include_cx=True
ne_cm3 – float or array Electron density in units of cm^-3
Te_eV – float or array Electron temperature in units of eV. If left to None, the Te grid given in the atomic data is used.
maxTe – float Maximum temperature of interest; only used if Te is left to None.
include_cx – bool If True, obtain charge exchange terms as well.

Returns

array (n_Te): log10 Te grid on which atomic rates are given
logS, logR (,logcx): arrays (n_ne,n_Te): atomic rates for effective ionization, radiative+dielectronic recombination (+ charge exchange, if requested). After exponentiation, all terms will be in units of s^-1.

Return type

logTe

aurora.atomic.get_frac_abundances(atom_data, ne_cm3, Te_eV=None, n0_by_ne=1e-05, include_cx=False, ne_tau=inf, plot=True, ax=None, rho=None, rho_lbl=None, ls='-')[source]¶

Calculate fractional abundances from ionization and recombination equilibrium. If include_cx=True, radiative recombination and thermal charge exchange are summed.

This method can work with ne,Te and n0_by_ne arrays of arbitrary dimension, but plotting is only supported in 1D (defaults to flattened arrays).

Parameters

atom_data – dictionary of atomic ADAS files (only acd, scd are required; ccd is necessary only if include_cx=True
ne_cm3 – float or array Electron density in units of cm^-3
Te_eV – float or array, optional Electron temperature in units of eV. If left to None, the Te grid given in the atomic data is used.
n0_by_ne – float or array, optional Ratio of background neutral hydrogen to electron density, used if include_cx=True.
include_cx – bool If True, charge exchange with background thermal neutrals is included.
ne_tau – float, opt Value of electron density in $m^{-3}\cdot s$ :math:` imes` particle residence time. This is a scalar value that can be used to model the effect of transport on ionization equilibrium. Setting ne_tau=np.inf (default) corresponds to no effect from transport.
plot – bool, optional Show fractional abundances as a function of ne,Te profiles parameterization.
ax – matplotlib.pyplot Axes instance Axes on which to plot if plot=True. If False, it creates new axes
rho – list or array, optional Vector of radial coordinates on which ne,Te (and possibly n0_by_ne) are given. This is only used for plotting, if given.
rho_lbl – str, optional Label to be used for rho. If left to None, defaults to a general “rho”.
ls – str, optional Line style for plots. Continuous lines are used by default.

Returns

array: log10 of electron temperatures as a function of which the fractional abundances and rate coefficients are given.
fzarray, (space,nZ): Fractional abundances across the same grid used by the input ne,Te values.
rate_coeffsarray, (space, nZ): Rate coefficients in units of [s^-1].

Return type

logTe

aurora.atomic.gff_mean(Z, Te)[source]¶: Total free-free gaunt factor yielding the total radiated bremsstrahlung power when multiplying with the result for gff=1. Data originally from Karzas & Latter, extracted from STRAHL’s atomic_data.f.

aurora.atomic.impurity_brems(nz, ne, Te)[source]¶

Approximate bremsstrahlung in units of $mW/nm/sr/m^3 \cdot cm^3$, or equivalently $W/m^3$ for full spherical emission.

Note that this may not be very useful, since this contribution is already included in the continuum radiation component in ADAS files. The bremsstrahlung estimate in ADAS continuum radiation files is more accurate than the one give by this function, which uses a simpler interpolation of the Gaunt factor with weak ne-dependence. Use with care!

Parameters

nz – array (time,nZ,space) Densities for each charge state [$cm^{-3}$]
ne – array (time,space) Electron density [$cm^{-3}$]
Te – array (time,space) Electron temperature [$cm^{-3}$]

Returns

array (time,nZ,space): Bremsstrahlung for each charge state

Return type

brems

aurora.atomic.interp_atom_prof(atom_table, x_prof, y_prof, log_val=False, x_multiply=True)[source]¶

Fast interpolate atomic data in atom_table onto the x_prof and y_prof profiles. This function assume that x_prof, y_prof, x,y, table are all base-10 logarithms, and x_prof, y_prof are equally spaced.

Parameters

atom_table – list List with x,y, table = atom_table, containing atomic data from one of the ADAS files.
x_prof – array (nt,nr) Spatio-temporal profiles of the first coordinate of the ADAS file table (usually electron density in cm^-3)
y_prof – array (nt,nr) Spatio-temporal profiles of the second coordinate of the ADAS file table (usually electron temperature in eV)
log_val – bool If True, return natural logarithm of the data
x_multiply – bool If True, multiply output by 10**x_prof.

Returns

array (nt,nion,nr): Interpolated atomic data on time,charge state and spatial grid that correspond to the ion of interest and the spatiotemporal grids of x_prof and y_prof.

Return type

interp_vals

aurora.atomic.null_space(A)[source]¶: Find null space of matrix A

aurora.atomic.plot_norm_ion_freq(S_z, q_prof, R_prof, imp_A, Ti_prof, nz_profs=None, rhop=None, plot=True, eps_prof=None)[source]¶

Compare effective ionization rate for each charge state with the characteristic transit time that a non-trapped and trapped impurity ion takes to travel a parallel distance L = q R.

If the normalized ionization rate is less than 1, then flux surface averaging of background asymmetries (e.g. from edge or beam neutrals) can be considered in a “flux-surface-averaged” sense; otherwise, local effects (i.e. not flux-surface-averaged) may be too important to ignore.

This function is inspired by Dux et al. NF 2020. Note that in this paper the ionization rate averaged over all charge state densities is considered. This function avoids the averaging over charge states, unless these are provided as an input.

Parameters

S_z – array (r,cs) [s^-1] Effective ionization rates for each charge state as a function of radius. Note that, for convenience within aurora, cs includes the neutral stage.
q_prof – array (r,) Radial profile of safety factor
R_prof – array (r,) or float [m] Radial profile of major radius, either given as an average of HFS and LFS, or also simply as a scalar (major radius on axis)
imp_A – float [amu] Atomic mass number, i.e. number of protons + neutrons (e.g. 2 for D)
Ti_prof – array (r,) Radial profile of ion temperature [eV]
nz_profs – array (r,cs), optional Radial profile for each charge state. If provided, calculate average normalized ionization rate over all charge states.
rhop – array (r,), optional Sqrt of poloidal flux radial grid. This is used only for (optional) plotting.
plot – bool, optional If True, plot results.
eps_prof – array (r,), optional Radial profile of inverse aspect ratio, i.e. r/R, only used if plotting is requested.

Returns

array (r,cs) or (r,): Normalized ionization rate. If nz_profs is given as an input, this is an average over all charge state; otherwise, it is given for each charge state.

Return type

nu_ioniz_star

aurora.atomic.plot_relax_time(logTe, rate_coeff, ax=None)[source]¶

Plot relaxation time of the ionization equilibrium corresponding to the inverse of the given rate coefficients.

Parameters

logTe – array (nr,) log-10 of Te [eV], on an arbitrary grid (same as other arguments, but not necessarily radial)
rate_coeff – array (nr,) Rate coefficients from ionization balance. See get_frac_abundances() to obtain these via the “compute_rates” argument. N.B.: these rate coefficients will depend also on electron density, which does affect relaxation times.
ax – matplotlib axes instance, optional If provided, plot relaxation times on these axes.

aurora.adas_files module¶

Functions to provide default ADAS files for Aurora modelling, including capabilities to fetch these files remotely from the OPEN-ADAS website.

aurora.adas_files.adas_files_dict()[source]¶

Selections for ADAS files for Aurora runs and radiation calculations. This function can be called to fetch a set of default files, which can then be modified (e.g. to use a new file for a specific SXR filter) before running a calculation.

Returns

dict: Dictionary with keys equal to the atomic symbols of many of the most common ions of interest in fusion research. For each ion, a sub-dictionary contains recommended file names for the relevant ADAS data types. Not all files types are available for all ions. Files types are usually a subset of ‘acd’,’scd’,’prb’,’plt’,’ccd’,’prc’,’pls’,’prs’,’fis’,’brs’,’pbs’,prc’ Refer to get_adas_file_types() for a description of the meaning of each file.

Return type

files

aurora.adas_files.fetch_adf11_file(filename)[source]¶

Download ADF11 file from the OPEN-ADAS website and store it in the ‘adas_data/adf11’ directory.

Parameters: filename – str Name of ADF11 file to be downloaded, e.g. ‘plt89_ar.dat’.

aurora.adas_files.fetch_adf15_file(filename)[source]¶

Download ADF15 file from the OPEN-ADAS website and store it in the ‘adas_data/adf15’ directory.

Parameters: filename – str Name of ADF15 file to be downloaded, e.g. ‘pec96#c_pju#c2.dat’.

aurora.adas_files.get_adas_file_loc(filename, filetype='adf11')[source]¶

Find location of requested atomic data file for the indicated ion. The search proceeds with the following attempts, in this order:

If the file is available in Aurora/adas_data/filetype, with filetype given by the user, always use this data.
If the environmental variable “AURORA_ADAS_DIR” is defined, attempt to find the file there and copy it to Aurora/adas_data/filetype.
Attempt to fetch the file remotely via open.adas.ac.uk and save it in Aurora/adas_data/filetype/.

Parameters

filename – str Name of the ADAS file of interest, e.g. ‘plt89_ar.dat’
filetype – str ADAS file type. Currently allowed: ‘adf11’ or ‘adf15’

Returns

str: Full path to the requested file.

Return type

file_loc

aurora.radiation module¶

aurora.radiation.adf04_files()[source]¶: Collection of trust-worthy ADAS ADF04 files. This function will be moved and expanded in ColRadPy in the near future.

aurora.radiation.compute_rad(imp, nz, ne, Te, n0=None, Ti=None, ni=None, adas_files_sub={}, prad_flag=False, sxr_flag=False, thermal_cx_rad_flag=False, spectral_brem_flag=False)[source]¶

Calculate radiation terms corresponding to a simulation result. The nz,ne,Te,n0,Ti,ni arrays are normally assumed to be given as a function of (time,nZ,space), but time and space may be substituted by other coordinates (e.g. R,Z)

Result can be conveniently plotted with a time-slider using, for example

aurora.slider_plot(rhop,time, res['line_rad'].transpose(1,2,0)/1e6,
    xlabel=r'$\rho_p$', ylabel='time [s]',
    zlabel=r'$P_{rad}$ [$MW$]',
    plot_sum=True,
    labels=[f'Ca$^{{{str(i)}}}$' for i in np.arange(res['line_rad'].shape[1])])

All radiation outputs are given in $W cm^{-3}$, consistently with units of $cm^{-3}$ given for inputs.

Parameters

imp – str Impurity symbol, e.g. Ca, F, W
nz – array (time, nZ, space) [$cm^{-3}$] Dictionary with impurity density result, as given by run_aurora() method.
ne – array (time,space) [$cm^{-3}$] Electron density on the output grids.
Te – array (time,space) [eV] Electron temperature on the output grids.

Keyword Arguments

n0 – array(time,space), optional [$cm^{-3}$] Background neutral density (assumed of hydrogen-isotopes). This is only used if thermal_cx_rad_flag=True.
Ti – array (time,space) [eV] Main ion temperature (assumed of hydrogen-isotopes). This is only used if thermal_cx_rad_flag=True. If not set, Ti is taken equal to Te.
adas_files_sub – dict Dictionary containing ADAS file names for radiation calculations, possibly including keys “plt”,”prb”,”prc”,”pls”,”prs”,”pbs”,”brs” Any file names that are needed and not provided will be searched in the adas_files_dict() dictionary.
prad_flag – bool, optional If True, total radiation is computed (for each charge state and their sum)
sxr_flag – bool, optional If True, soft x-ray radiation is computed (for the given ‘pls’,’prs’ ADAS files)
thermal_cx_rad_flag – bool, optional If True, thermal charge exchange radiation is computed.
spectral_brem_flag – bool, optional If True, spectral bremstrahlung is computed (based on available ‘brs’ ADAS file)

Returns

dict: Dictionary containing radiation terms, depending on the activated flags. The structure of the “res” dictionary is as follows.

If prad_flag=True,

res[‘line_rad’]array (nt,nZ,nr)- from ADAS “plt” files: Excitation-driven line radiation for each impurity charge state.
res[‘cont_rad’]array (nt,nZ,nr)- from ADAS “prb” files: Continuum and line power driven by recombination and bremsstrahlung for impurity ions.
res[‘brems’]array (nt,nr)- analytic formula.: Bremsstrahlung produced by electron scarrering at fully ionized impurity This is only an approximate calculation and is more accurately accounted for in the ‘cont_rad’ component.
res[‘thermal_cx_cont_rad’]array (nt,nZ,nr)- from ADAS “prc” files: Radiation deriving from charge transfer from thermal neutral hydrogen to impurity ions. Returned only if thermal_cx_rad_flag=True.
res[‘tot’]array (nt,nZ,nr): Total unfilted radiation, summed over all charge states, given by the sum of all known radiation components.

If sxr_flag=True,

res[‘sxr_line_rad’]array (nt,nZ,nr)- from ADAS “pls” files: Excitation-driven line radiation for each impurity charge state in the SXR range.
res[‘sxr_cont_rad’]array (nt,nZ,nr)- from ADAS “prs” files: Continuum and line power driven by recombination and bremsstrahlung for impurity ions in the SXR range.
res[‘sxr_brems’]array (nt,nZ,nr)- from ADAS “pbs” files: Bremsstrahlung produced by electron scarrering at fully ionized impurity in the SXR range.
res[‘sxr_tot’]array (nt,nZ,nr): Total radiation in the SXR range, summed over all charge states, given by the sum of all known radiation components in the SXR range.

If spectral_brem_flag,

res[‘spectral_brems’]array (nt,nZ,nr) – from ADAS “brs” files: Bremsstrahlung at a specific wavelength, depending on provided “brs” file.

Return type

res

aurora.radiation.get_colradpy_pec_prof(ion, cs, rhop, ne_cm3, Te_eV, lam_nm=1.8705, lam_width_nm=0.002, meta_idxs=[0], adf04_repo='/home/docs/adf04_files/ca/ca_adf04_adas/', pec_threshold=1e-20, phot2energy=True, plot=True)[source]¶

Compute radial profile for Photon Emissivity Coefficients (PEC) for lines within the chosen wavelength range using the ColRadPy package. This is an alternative to the option of using the read_adf15() function to read PEC data from an ADAS ADF-15 file and interpolate results on ne,Te grids.

Parameters

ion – str Ion atomic symbol
cs – str Charge state, given in format like ‘Ca18+’
rhop – array (nr,) Srt of normalized poloidal flux radial array
ne_cm3 – array (nr,) Electron density in $cm^{-3}$ units
Te_eV – array (nr,) Electron temperature in eV units
lam_nm – float Center of the wavelength region of interest [nm]
lam_width_nm – float Width of the wavelength region of interest [nm]
meta_idxs – list of integers List of levels in ADF04 file to be treated as metastable states.
adf04_repo – str Location where ADF04 file from :py:method:adf04_files() should be fetched.
prec_threshold – float Minimum value of PECs to be considered, in $photons \cdot cm^3/s$
phot2energy – bool If True, results are converted from $photons \cdot cm^3/s$ to $W.cm^3$
plot – bool If True, plot lines profiles and total

Returns

array (nr,): Radial profile of PEC intensity, in units of $photons \cdot cm^3/s$ (if phot2energy=False) or $W \cdot cm^3$ depending (if phot2energy=True).

Return type

pec_tot_prof

aurora.radiation.get_main_ion_dens(ne_cm3, ions, rhop_plot=None)[source]¶

Estimate the main ion density via quasi-neutrality. This requires subtracting from ne the impurity charge state density times Z for each charge state of every impurity present in the plasma in significant amounts.

Parameters

ne_cm3 – array (time,space) Electron density in $cm^{-3}$
ions – dict Dictionary with keys corresponding to the atomic symbol of each impurity under consideration. The values in ions[key] should correspond to the charge state densities for the selected impurity ion in $cm^{-3}$, with shape (time,nZ,space).
rhop_plot – array (space), optional rhop radial grid on which densities are given. If provided, plot densities of all species at the last time slice over this radial grid.

Returns

array (time,space): Estimated main ion density in $cm^{-3}$.

Return type

ni_cm3

aurora.radiation.plot_radiation_profs(imp, nz_prof, logne_prof, logTe_prof, xvar_prof, xvar_label='', atom_data=None)[source]¶

Compute profiles of predicted radiation, both SXR-filtered and unfiltered. This function offers a simplified interface to radiation calculation with respect to compute_rad(), which is more complete.

This function can be used to plot radial profiles (setting xvar_prof to a radial grid) or profiles as a function of any variable on which the logne_prof and logTe_prof may depend.

The variable “nz_prof” may be a full description of impurity charge state densities (e.g. the output of aurora), or profiles of fractional abundances from ionization equilibrium.

Parameters

imp – str, optional Impurity ion atomic symbol.
nz_prof – array (TODO for docs: check dimensions) Impurity charge state densities
logne_prof – array (TODO for docs: check dimensions) Electron density profiles in $cm^{-3}$
logTe_prof – array (TODO for docs: check dimensions) Electron temperature profiles in eV
xvar_prof – array (TODO for docs: check dimensions) Profiles of a variable of interest, on the same grid as kinetic profiles.
xvar_label – str, optional Label for x-axis.
atom_data – dict, optional Dictionary containing atomic data as output by get_atom_data() for the atomic processes of interest. “prs”,”pls”,”plt” and “prb” are required by this function. If not provided, this function loads these files internally.

Returns

array (TODO for docs: check dimensions): SXR line radiation.
prsarray (TODO for docs: check dimensions): SXR continuum radiation.
plttarray (TODO for docs: check dimensions): Unfiltered line radiation.
prbarray (TODO for docs: check dimensions): Unfiltered continuum radiation.

Return type

pls

aurora.radiation.radiation_model(imp, rhop, ne_cm3, Te_eV, vol, adas_files_sub={}, n0_cm3=None, Ti_eV=None, nz_cm3=None, frac=None, plot=False)[source]¶

Model radiation from a fixed-impurity-fraction model or from detailed impurity density profiles for the chosen ion. This method acts as a wrapper for :py:method:compute_rad(), calculating radiation terms over the radius and integrated over the plasma cross section.

Parameters

imp – str (nr,) Impurity ion symbol, e.g. W
rhop – array (nr,) Sqrt of normalized poloidal flux array from the axis outwards
ne_cm3 – array (nr,) Electron density in $cm^{-3}$ units.
Te_eV – array (nr,) Electron temperature in eV
vol – array (nr,) Volume of each flux surface in $m^3$. Note the units! We use $m^3$ here rather than $cm^3$ because it is more common to work with $m^3$ for flux surface volumes of fusion devices.

Keyword Arguments

adas_files_sub – dict Dictionary containing ADAS file names for forward modeling and/or radiation calculations. Possibly useful keys include “scd”,”acd”,”ccd”,”plt”,”prb”,”prc”,”pls”,”prs”,”pbs”,”brs” Any file names that are needed and not provided will be searched in the adas_files_dict() dictionary.
n0_cm3 – array (nr,), optional Background ion density (H,D or T). If provided, charge exchange (CX) recombination is included in the calculation of charge state fractional abundances.
Ti_eV – array (nr,), optional Background ion density (H,D or T). This is only used if CX recombination is requested, i.e. if n0_cm3 is not None. If not given, Ti is set equal to Te.
nz_cm3 – array (nr,nz), optional Impurity charge state densities in $cm^{-3}$ units. Fractional abundancies can alternatively be specified via the :param:frac parameter for a constant-fraction impurity model across the radius. If provided, nz_cm3 is used.
frac – float, optional Fractional abundance, with respect to ne, of the chosen impurity. The same fraction is assumed across the radial profile. If left to None, nz_cm3 must be given.
plot – bool, optional If True, plot a number of diagnostic figures.

Returns

dict: Dictionary containing results of radiation model.

Return type

res

aurora.radiation.read_adf15(path, order=1, plot_lines=[], ax=None, Te_max=None, ne_max=None, plot_log=False, plot_3d=False, pec_plot_min=None, pec_plot_max=None)[source]¶

Read photon emissivity coefficients from an ADAS ADF15 file.

Returns a dictionary whose keys are the wavelengths of the lines in angstroms. The value is an interp2d instance that will evaluate the PEC at a desired density and temperature.

Parameters

path – str Path to adf15 file to read.
order – int, opt Parameter to control the order of interpolation.

Keyword Arguments

plot_lines – list List of lines whose PEC data should be displayed. Lines should be identified by their wavelengths. The list of available wavelengths in a given file can be retrieved by first running this function ones, checking dictionary keys, and then requesting a plot of one (or more) of them.
ax – matplotlib axes instance If not None, plot on this set of axes
plot_log – bool When plotting, set a log scale
plot_3d – bool Display PEC data as a 3D plot rather than a 2D one.
pec_plot_min – float Minimum value of PEC to visualize in a plot
pec_plot_max – float Maximum value of PEC to visualize in a plot
Te_max – float Maximum Te value to plot when len(plot_lines)>1
ne_max – float Maximum ne value to plot when len(plot_lines)>1

Returns

dict: Dictionary containing interpolation functions for each of the available lines of the indicated type (ionization or recombination). Each interpolation function takes as arguments the log-10 of ne and Te.

Return type

pec_dict

Minimal Working Example (MWE):

filename = 'pec96#h_pju#h0.dat' # for D Ly-alpha

# fetch file automatically, locally, from AURORA_ADAS_DIR, or directly from the web:
path = aurora.get_adas_file_loc(filename, filetype='adf15')

# plot Lyman-alpha line at 1215.2 A. See available lines with pec_dict.keys() after calling without plot_lines argument
pec_dict = aurora.read_adf15(path, plot_lines=[1215.2])

Another example, this time also with charge exchange::

filename = 'pec96#c_pju#c2.dat'
path = aurora.get_adas_file_loc(filename, filetype='adf15')
pec_dict = aurora.read_adf15(path, plot_lines=[361.7])

Metastable-resolved files will be automatically identified and parsed accordingly, e.g.

filename = ‘pec96#he_pjr#he0.dat’ path = aurora.get_adas_file_loc(filename, filetype=’adf15’) pec_dict = aurora.read_adf15(path, plot_lines=[584.4])

This function should work with PEC files produced via adas810 or adas218.

aurora.grids_utils module¶

Methods to create radial and time grids for aurora simulations.

aurora.grids_utils.create_aurora_time_grid(timing, plot=False)[source]¶

Create time grid for simulations using a Fortran routine for definitions. The same functionality is offered by create_time_grid(), which however is written in Python. This method is legacy code; it is recommended to use the other.

Parameters

timing – dict Dictionary containing timing[‘times’],timing[‘dt_start’],timing[‘steps_per_cycle’],timing[‘dt_increase’] which define the start times to change dt values at, the dt values to start with, the number of time steps before increasing the dt by dt_increase. The last value in each of these arrays is used for sawteeth, whenever these are modelled, or else are ignored. This is the same time grid definition as used in STRAHL.
plot – bool, optional If True, display the created time grid.

Returns

array: Computational time grid corresponding to timing input.
savearray: Array of zeros and ones, where ones indicate that the time step will be stored in memory in aurora simulations. Points corresponding to zeros will not be returned to spare memory.

Return type

time

aurora.grids_utils.create_radial_grid(namelist, plot=False)[source]¶

Create radial grid for Aurora based on K, dr_0, dr_1, rvol_lcfs and bound_sep parameters. The lim_sep parameters is additionally used if plotting is requested.

Radial mesh points are set to be equidistant in the coordinate $\rho$, with

\[\rho = \frac{r}{\Delta r_{centre}} + \frac{r_{edge}}{k+1} \left(\frac{1}{\Delta r_{edge}}- \frac{1}{\Delta r_{centre}} \right) \left(\frac{r}{r_{edge}} \right)^{k+1}\]

The corresponding radial step size is

\[\Delta r = \left[\frac{1}{\Delta r_{centre}} + \left(\frac{1}{\Delta r_{edge}} - \frac{1}{\Delta r_{centre}} \right) \left(\frac{r}{r_{edge}}\right)^k \right]^{-1}\]

See the STRAHL manual for details.

Parameters

namelist – dict Dictionary containing Aurora namelist. This function uses the K, dr_0, dr_1, rvol_lcfs and bound_sep parameters. Additionally, lim_sep is used if plotting is requested.
plot – bool, optional If True, plot the radial grid spacing vs. radial location.

Returns

array: Volume-normalized grid used for Aurora simulations.
proarray: Normalized first derivatives of the radial grid, defined as pro = (drho/dr)/(2 d_rho) = rho’/(2 d_rho)
qprarray: Normalized second derivatives of the radial grid, defined as qpr = (d^2 rho/dr^2)/(2 d_rho) = rho’‘/(2 d_rho)
prox_paramfloat: Grid parameter used for perpendicular loss rate at the last radial grid point.

Return type

rvol_grid

aurora.grids_utils.create_time_grid(timing=None, plot=False)[source]¶

Create time grid for simulations using the Fortran implementation of the time grid generator.

Parameters

timing – dict Dictionary containing timing elements: ‘times’, ‘dt_start’, ‘steps_per_cycle’,’dt_increase’ As in STRAHL, the last element in each of these arrays refers to sawtooth events.
plot – bool If True, plot time grid.

Returns

array: Computational time grid corresponding to :param:timing input.
savearray: Array of zeros and ones, where ones indicate that the time step will be stored in memory in Aurora simulations. Points corresponding to zeros will not be returned to spare memory.

Return type

time

aurora.grids_utils.create_time_grid_new(timing, verbose=False, plot=False)[source]¶

Define time base for Aurora based on user inputs This function reproduces the functionality of STRAHL’s time_steps.f Refer to the STRAHL manual for definitions of the time grid

Parameters

n – int Number of elements in time definition arrays
t – array Time vector of the time base changes
dtstart – array dt value at the start of a cycle
itz – array cycle length, i.e. number of time steps before increasing dt
tinc – factor by which time steps should be increasing within a cycle
verbose – bool If Trueprint to terminal a few extra info

Returns

array: Times in the time base [s]
i_savearray: Array of 0,1 values indicating at which times internal arrays should be stored/returned.

Return type

t_vals

~~~~~~~~~~~ THIS ISN’T FUNCTIONAL YET! ~~~~~~~~~~~~

aurora.grids_utils.estimate_boundary_distance(shot, device, time_ms)[source]¶

Obtain a simple estimate for the distance between the LCFS and the wall boundary. This requires access to the A_EQDSK on the EFIT01 tree on MDS+. Users who may find that this call does not work for their device may try to adapt the OMFITmdsValue TDI string.

Parameters

shot – int Discharge/experiment number
device – str Name of device, e.g. ‘C-Mod’, ‘DIII-D’, etc.
time_ms – int or float Time at which results for the outer gap should be taken.

Returns

float: Estimate for the distance between the wall boundary and the separatrix [cm]
lim_sepfloat: Estimate for the distance between the limiter and the separatrix [cm]. This is (quite arbitrarily) taken to be 2/3 of the bound_sep distance.

Return type

bound_sep

aurora.grids_utils.estimate_clen(geqdsk)[source]¶

Estimate average connection length in the open SOL and in the limiter shadow NB: these are just rough numbers!

Parameters

geqdsk – dict EFIT g-EQDSK as processed by the omfit_eqdsk package.

Returns

float: Estimate of the connection length to the divertor
clen_limiterfloat: Estimate of the connection length to the limiter

Return type

clen_divertor

aurora.grids_utils.get_HFS_LFS(geqdsk, rho_pol=None)[source]¶

Get high-field-side (HFS) and low-field-side (LFS) major radii from the g-EQDSK data. This is useful to define the rvol grid outside of the LCFS. See the get_rhopol_rV_mapping() for an application.

Parameters

geqdsk – dict Dictionary containing the g-EQDSK file as processed by the omfit_eqdsk package.
rho_pol – array, optional Array corresponding to a grid in sqrt of normalized poloidal flux for which a corresponding rvol grid should be found. If left to None, an arbitrary grid will be created internally.

Returns

array: Major radius [m] on the HFS
Rlfsarray: Major radius [m] on the LFS

Return type

Rhfs

aurora.grids_utils.get_rhopol_rvol_mapping(geqdsk, rho_pol=None)[source]¶

Compute arrays allowing 1-to-1 mapping of rho_pol and rvol, both inside and outside the LCFS.

rvol is defined as $\sqrt{V/(2 \pi^2 R_{axis}}$ inside the LCFS. Outside of it, we artificially expand the LCFS to fit true equilibrium at the midplane based on the rho_pol grid (sqrt of normalized poloidal flux).

Method:

\[\begin{split}r(\rho,\theta) = r_0(\rho) + (r_{lcfs}(\theta) - r_{0,lcfs}) \times \mathcal{f} \\ z(\rho,\theta) = z_0 + (z_{lcfs}(\theta) - z_0 ) \times \mathcal{f} \\ \mathcal{f} = \frac{ r(\rho,\theta=0) - r(\rho,\theta=180)}{r_{lcfs}(\theta=0)- r_{lcfs}(\theta=180)} \\ r_{0,lcfs} = \frac{1}{2} (r_{lcfs}(\theta=0)+ r_{lcfs}(\theta=180)) \\ r_0(\rho) = \frac{1}{2} (r(\rho,\theta=0) + r(\rho,\theta=180))\end{split}\]

The mapping between rho_pol and rvol allows one to interpolate inputs on a rho_pol grid onto the rvol grid (in cm) used internally by the code.

Parameters

geqdsk – dict Dictionary containing the g-EQDSK file as processed by the omfit_eqdsk package.
rho_pol – array, optional Array corresponding to a grid in sqrt of normalized poloidal flux for which a corresponding rvol grid should be found. If left to None, an arbitrary grid will be created internally.

Returns

array: Sqrt of normalized poloidal flux grid
rvolarray: Mapping of rho_pol to a radial grid defined in terms of normalized flux surface volume.

Return type

rho_pol

aurora.coords module¶

aurora.coords.rV_vol_average(quant, r_V)[source]¶

Calculate a volume average of the given radially-dependent quantity on a r_V grid.

This function makes useof the fact that the r_V radial coordinate, defined as r_V = sqrt{ V / (2 pi^2 R_{axis} }, maps shaped volumes onto a circular geometry, making volume averaging a trivial operation via langle Q

angle = Sigma_i Q(r_i) 2 pi Delta r_V

where $Delta r_V$ is the spacing between radial points in r_V.

Note that if the input r_V coordinate is extended outside the LCFS, this function will return the effective volume average also in the SOL, since it is agnostic to the presence of the LCFS.

Args:

quantarray, (space, …): quantity that one wishes to volume-average. The first dimension must correspond to r_V, but other dimensions may be exist afterwards.
r_Varray, (space,): Radial r_V coordinate in cm units.

Returns:

quant_vol_avgarray, (space, …): Volume average of the quantity given as an input, in the same units as in the input

aurora.coords.rad_coord_transform(x, name_in, name_out, geqdsk)[source]¶

Transform from one radial coordinate to another. Note that this coordinate conversion is only strictly valid inside of the LCFS.

Parameters

x – array input x coordinate
name_in – str input x coordinate name (‘rhon’,’r_V’,’rhop’,’rhov’,’Rmid’,’rmid’,’roa’)
name_out – str input x coordinate (‘rhon’, ‘r_V’, ‘rhop’,’rhov’,’Rmid’,’rmid’,’roa’)
geqdsk – dict gEQDSK dictionary, as obtained from the omfit-eqdsk package.

Returns

Conversion of x for the requested radial grid coordinate.

aurora.coords.vol_average(quant, rhop, method='omfit', geqdsk=None, device=None, shot=None, time=None, return_geqdsk=False)[source]¶

Calculate the volume average of the given radially-dependent quantity on a rhop grid.

Parameters

quant – array, (space, …) quantity that one wishes to volume-average. The first dimension must correspond to space, but other dimensions may be exist afterwards.
rhop – array, (space,) Radial rhop coordinate in cm units.
method – {‘omfit’,’fs’} Method to evaluate the volume average. The two options correspond to the way to compute volume averages via the OMFIT fluxSurfaces classes and via a simpler cumulative sum in r_V coordinates. The methods only slightly differ in their results. Note that ‘omfit’ will fail if rhop extends beyond the LCFS, while method ‘fs’ can estimate volume averages also into the SOL. Default is method=’omfit’.
geqdsk – output of the omfit_eqdsk.OMFITgeqdsk class, postprocessing the EFIT geqdsk file containing the magnetic geometry. If this is left to None, the function internally tries to fetch it using MDS+ and omfit_eqdsk. In this case, device, shot and time to fetch the equilibrium are required.
device – str Device name. Note that routines for this device must be implemented in omfit_eqdsk for this to work.
shot – int Shot number of the above device, e.g. 1101014019 for C-Mod.
time – float Time at which equilibrium should be fetched in units of ms.
return_geqdsk – bool If True, omfit_eqdsk dictionary is also returned

Returns

array, (space, …): Volume average of the quantity given as an input, in the same units as in the input. If extrapolation beyond the range available from EFIT volume averages over a shorter section of the radial grid will be attempted. This does not affect volume averages within the LCFS.
geqdskdict: Only returned if return_geqdsk=True.

Return type

quant_vol_avg

aurora.source_utils module¶

Methods related to impurity source functions.

sciortino, 2020

aurora.source_utils.get_radial_source(namelist, rvol_grid, pro_grid, S_rates, Ti_eV=None)[source]¶

Obtain spatial dependence of source function.

If namelist[‘source_width_in’]==0 and namelist[‘source_width_out’]==0, the source radial profile is defined as an exponential decay due to ionization of neutrals. This requires S_rates, the ionization rate of neutral impurities, to be given with S_rates.shape=(len(rvol_grid),len(time_grid))

If namelist[‘imp_source_energy_eV’]<0, the neutrals speed is taken as the thermal speed based on Ti_eV, otherwise the value corresponding to namelist[‘imp_source_energy_eV’] is used.

Parameters

namelist – dict Aurora namelist. Only elements referring to the spatial distribution and energy of source atoms are accessed.
rvol_grid – array (nr,) Radial grid in volume-normalized coordinates [cm]
pro_grid – array (nr,) Normalized first derivatives of the radial grid in volume-normalized coordinates.
S_rates – array (nr,nt) Ionization rate of neutral impurity over space and time.

Keyword Arguments

Ti_eV – array, optional (nr,nt) Background ion temperature, only used if source_width_in=source_width_out=0.0 and imp_source_energy_eV<=0, in which case the source impurity neutrals are taken to have energy equal to the local Ti [eV].

Returns

array (nr,nt): Radial profile of the impurity neutral source for each time step.

Return type

source_rad_prof

aurora.source_utils.get_source_time_history(namelist, Raxis_cm, time)[source]¶

Load source time history based on current state of the namelist.

There are 4 options to describe the time-dependence of the source:

(1) namelist[‘source_type’] == ‘file’: in this case, a simply formatted source file, with one time point and corresponding and source amplitude on each line, is read in. This can describe an arbitrary time dependence, e.g. as measured from an experimental diagnostic.

(2) namelist[‘source_type’] == ‘const’: in this case, a constant source (e.g. a gas puff) is simulated. It is recommended to run the simulation for >100ms in order to see self-similar charge state profiles in time.

(3) namelist[‘source_type’] == ‘step’: this allows the creation of a source that suddenly appears and suddenly stops, i.e. a rectangular “step”. The duration of this step is given by namelist[‘step_source_duration’]. Multiple step times can be given as a list in namelist[‘src_step_times’]; the amplitude of the source at each step is given in namelist[‘src_step_rates’]

(4) namelist[‘source_type’] == ‘synth_LBO’: this produces a model source from a LBO injection, given by a convolution of a gaussian and an exponential. The required parameters in this case are inside a namelist[‘LBO’] dictionary: namelist[‘LBO’][‘t_start’], namelist[‘LBO’][‘t_rise’], namelist[‘LBO’][‘t_fall’], namelist[‘LBO’][‘n_particles’]. The “n_particles” parameter corresponds to the amplitude of the source (the number of particles corresponding to the integral over the source function.

Parameters

namelist – dict Aurora namelist dictionary.
Raxis_cm – float Major radius at the magnetic axis [cm]. This is needed to normalize the source such that it is treated as toroidally symmetric – a necessary idealization for 1.5D simulations.
time – array (nt,), optional Time array the source should be returned on.

Returns

array (nt,): The source time history on the input time base.

Return type

source_time_history

aurora.source_utils.lbo_source_function(t_start, t_rise, t_fall, n_particles=1.0, time_vec=None)[source]¶

Model for the expected shape of the time-dependent source function, using a convolution of a gaussian and an exponential decay.

Parameters

t_start – float or array-like [ms] Injection time, beginning of source rise. If multiple values are given, they are used to create multiple source functions.
t_rise – float or array-like [ms] Time scale of source rise. Similarly to t_start for multiple values.
t_fall – float or array-like [ms] Time scale of source decay.Similarly to t_start for multiple values.
n_particles – float, opt Total number of particles in source. Similarly to t_start for multiple values. Defaults to 1.0.
time_vec – array-like Time vector on which to create source function. If left to None, use a linearly spaced time vector including the main features of the function.

Returns

array: Times for the source function of each given impurity
sourcearray: Time history of the synthetized source function.

Return type

time_vec

aurora.source_utils.read_source(filename)[source]¶

Read a STRAHL source file from {imp}flx{shot}.dat locally.

Parameters

filename – str Location of the file containing the STRAHL source file.

Returns

array of float, (n,): The timebase (in seconds).
sarray of float, (n,): The source function (#/s).

Return type

t

aurora.source_utils.write_source(t, s, shot, imp='Ca')[source]¶

Write a STRAHL source file.

This will overwrite any {imp}flx{shot}.dat locally.

Parameters

t – array of float, (n,) The timebase (in seconds).
s – array of float, (n,) The source function (in particles/s).
shot – int Shot number, only used for saving to a .dat file
imp – str, optional Impurity species atomic symbol

Returns

str: Content of the source file written to {imp}flx{shot}.dat

Return type

contents

aurora.plot_tools module¶

aurora.plot_tools.get_color_cycle()[source]¶

aurora.plot_tools.get_line_cycle()[source]¶

aurora.plot_tools.get_ls_cycle()[source]¶

aurora.plot_tools.slider_plot(x, y, z, xlabel='', ylabel='', zlabel='', labels=None, plot_sum=False, x_line=None, y_line=None, **kwargs)[source]¶

Make a plot to explore multidimensional data.

Parameters

x – array of float, (M,) The abscissa. (in aurora, often this may be rhop)
y – array of float, (N,) The variable to slide over. (in aurora, often this may be time)
z – array of float, (P, M, N) The variables to plot.
xlabel – str, optional The label for the abscissa.
ylabel – str, optional The label for the slider.
zlabel – str, optional The label for the ordinate.
labels – list of str with length P The labels for each curve in z.
plot_sum – bool, optional If True, will also plot the sum over all P cases. Default is False.
x_line – float, optional x coordinate at which a vertical line will be drawn.
y_line – float, optional y coordinate at which a horizontal line will be drawn.

aurora.default_nml module¶

Method to load default namelist. This should be complemented with additional info by each user.

sciortino, July 2020

aurora.default_nml.load_default_namelist()[source]¶: Load default namelist. Users should modify and complement this for a successful run.

aurora.interp module¶

This script contains a number of functions used for interpolation of kinetic profiles and D,V profiles in STRAHL. Refer to the STRAHL manual for details.

aurora.interp.exppol0(params, d, rLCFS, r)[source]¶

aurora.interp.exppol1(params, d, rLCFS, r)[source]¶

aurora.interp.funct(params, rLCFS, r)[source]¶

Function ‘funct’ in STRAHL manual

The “params” input is broken down into 6 arguments:: y0 is core offset y1 is edge offset y2 (>y0, >y1) sets the gaussian amplification p0 sets the width of the inner gaussian P1 sets the width of the outer gaussian p2 sets the location of the inner and outer peaks

aurora.interp.funct2(params, rLCFS, r)[source]¶: Function ‘funct2’ in STRAHL manual.

aurora.interp.interp(x, y, rLCFS, r)[source]¶: Function ‘interp’ used in STRAHL for D and V profiles.

aurora.interp.interp_quad(x, y, d, rLCFS, r)[source]¶: Function ‘interp’ used for kinetic profiles.

aurora.interp.interpa_quad(x, y, rLCFS, r)[source]¶: Function ‘interpa’ used for kinetic profiles

aurora.interp.ratfun(params, d, rLCFS, r)[source]¶

aurora.animate module¶

aurora.animate.animate_aurora(x, y, z, xlabel='', ylabel='', zlabel='', labels=None, plot_sum=False, uniform_y_spacing=True, save_filename=None)[source]¶

Produce animation of time- and radially-dependent results from aurora.

Parameters

x – array of float, (M,) The abscissa. (in aurora, often this may be rhop)
y – array of float, (N,) The variable to slide over. (in aurora, often this may be time)
z – array of float, (P, M, N) The variables to plot.
xlabel – str, optional The label for the abscissa.
ylabel – str, optional The label for the animated coordinate. This is expected in a format such that ylabel.format(y_val) will display a good moving label, e.g. ylabel=’t={:.4f} s’.
zlabel – str, optional The label for the ordinate.
labels – list of str with length P The labels for each curve in z.
plot_sum – bool, optional If True, will also plot the sum over all P cases. Default is False.
uniform_y_spacing – bool, optional If True, interpolate values in z onto a uniformly-spaced y grid
save_filename – str If a valid path/filename is provided, the animation will be saved here in mp4 format.

aurora.particle_conserv module¶

aurora.particle_conserv.check_particle_conserv(Raxis_cm, ds=None, filepath=None, linestyle='-', plot=True, axs=None)[source]¶

Check time evolution and particle conservation in Aurora or STRAHL output.

Parameters

Raxis_cm – float Major radius on axis [cm], used for volume integrals.
ds – xarray dataset, optional Dataset containing Aurora results, created using the xarray package. See check_conservation() for an illustration on how to use this.
filepath – str, optional If provided, load results from STRAHL output file and check particle particle conservation as for an Aurora run.
linestyle – str, optional matplotlib linestyle, default is ‘-‘ (continuous lines). Use this to overplot lines on the same plots using different linestyles, e.g. to check whether some aurora option makes particle conservation better or worse.
plot – bool, optional If True, plot time histories of particle densities in each simulation reservoir.
axs – 2-tuple or array array-like structure containing two matplotlib.Axes instances: the first one for the separate particle time variation in each reservoir, the second for the total particle-conservation check. This can be used to plot results from several aurora runs on the same axes.

Returns

dict: Dictionary containing time histories across all reservoirs, useful for the assessment of particle conservation.
axs2-tuple or array, only returned if plot=True: array-like structure containing two matplotlib.Axes instances, (ax1,ax2). See optional input argument.

Return type

out

aurora.particle_conserv.vol_int(Raxis_cm, ds, var, rhop_max=None)[source]¶

Perform a volume integral of an input variable. If the variable is f(t,x) then the result is f(t). If the variable is f(t,*,x) then the result is f(t,charge) when “*” represents charge, line index, etc…

Parameters

Raxis_cm – float Major radius on axis [cm]
ds – xarray dataset Dataset containing Aurora or STRAHL result
var – str Name of the variable in the strahl_result.cdf file
rhop_max – float Maximum normalized poloidal flux for integral. If not provided, integrate over the entire simulation grid.

Returns

array (nt,): Time history of volume integrated variable

Return type

var_volint

aurora.neutrals module¶

Aurora functionality for edge neutral modeling. The ehr5 file from DEGAS2 is used. See https://w3.pppl.gov/degas2/ for details.

aurora.neutrals.download_ehr5_file()[source]¶

Download the ehr5.dat file containing atomic data describing the multi-step ionization and recombination of hydrogen.

See https://w3.pppl.gov/degas2/ for details.

class aurora.neutrals.ehr5_file(filepath=None)[source]¶

Bases: object

Read ehr5.dat file from DEGAS2. Returns a dictionary containing

Ionization rate Seff in $cm^3 s^{-1}$
Recombination rate Reff in $cm^3 s^{-1}$
Neutral electron losses $E_{loss}^{(i)}$ in $erg s^{-1}$
Continuum electron losses $E_{loss}^{(ii)}$ in $erg s^{-1}$
Neutral “n=2 / n=1”, $N_2^{(i)}/N_1$
Continuum “n=2 / n=1”, :math:`N_2^{(ii)}/N_11
Neutral “n=3 / n=1”, $N_3^{(i)}/N_1$
Continuum “n=3 / n=1”, $N_3^{(ii)}/N_1$

… and similarly for n=4 to 9.

Refer to the DEGAS2 manual for details.

load()[source]¶

plot(field='Seff', fig=None, axes=None)[source]¶

aurora.neutrals.get_exc_state_ratio(m, N1, ni, ne, Te, rad_prof=None, rad_label='rmin [cm]', plot=False)[source]¶

Compute density of excited states in state m (m>1), given the density of ground state atoms.

This function is not l-resolved.

The function returns

\[N_m/N_1 = \left(\]

rac{N_m^i}{N_1} ight) N_m + left( rac{N_m^{ii}}{n_i} ight) n_i

where $N_m$ is the number of electrons in the excited state $m$, $N_1$ is the number in the ground state, and $n_i$ is the density of ions that could recombine. $i$ and $ii$ indicate terms corresponding to coupling to the ground state and to the continuum, respectively.

Ref.: DEGAS2 manual.

Args:

mint
Principal quantum number of excited state of interest. 2<m<10

N1float, list or 1D-array [$cm^-3$]
Density of ions in the ground state. This must have the same shape as ni!

nifloat, list or 1D-array [$cm^-3$]
Density of ions corresponding to the atom under consideration. This must have the same shape as N1!

nefloat, list or 1D-array [$cm^-3$]
Electron density to evaluate atomic rates at.

Tefloat, list or 1D-array [$eV$]
Electron temperature to evaluate atomic rates at.

Keyword Args:

rad_proflist, 1D array or None
If None, excited state densities are evaluated at all the combinations of ne,Te and zip(Ni,ni). If a 1D array (same length as ne,Te,ni and N1), then this is taken to be a radial coordinate for radial profiles of ne,Te,ni and N1.

rad_labelstr
When rad_prof is not None, this is the label for the radial coordinate.

plotbool
Display the excited state ratio

Returns:

Nmarray of shape [len(ni)=len(N1),len(ne),len(Te)]
Density of electrons in excited state n [$cm^-3$]

aurora.neutrals.plot_exc_ratios(n_list=[2, 3, 4, 5, 6, 7, 8, 9], ne=10000000000000.0, ni=10000000000000.0, Te=50, N1=1000000000000.0, ax=None, ls='-', c='r', label=None)[source]¶

Plot $N_i/N_1$, the ratio of hydrogen neutral density in the excited state i and the ground state, for given electron density and temperature.

Parameters

n_list – list of integers List of excited states (principal quantum numbers) to consider.
ne – float Electron density in $cm^-3$.
ni – float Ionized hydrogen density [$cm^-3$]. This may be set equal to ne for a pure plasma.
Te – float Electron temperature in $eV$.
N1 – float Density of ground state hydrogen [$cm^-3$]. This is needed because the excited state fractions depend on the balance of excitation from the ground state and coupling to the continuum.

Keyword Arguments

ax – matplotlib.axes instance, optional Axes instance on which results should be plotted.
ls – str Line style to use
c – str or other matplotlib color specification Color to use in plots
label – str Label to use in scatter plot.

Returns

list of arrays: List of arrays for each of the n-levels requested, each containing excited state densities at the chosen densities and temperatures for the given ground state density.

Return type

Ns

aurora.nbi_neutrals module¶

Methods for neutral beam analysis, particularly in relation to impurity transport studies. These script collects functions that should be device-agnostic.

aurora.nbi_neutrals.beam_grid(uvw_src, axis, max_radius=255.0)[source]¶

Method to obtain the 3D orientation of a beam with respect to the device. The uvw_src and (normalized) axis arrays may be obtained from the d3d_beams method of fidasim_lib.py in the FIDASIM module in OMFIT.

This is inspired by beam_grid in fidasim_lib.py of the FIDASIM module (S. Haskey) in OMFIT.

aurora.nbi_neutrals.bt_rate_maxwell_average(sigma_fun, Ti, E_beam, m_bckg, m_beam, n_level)[source]¶

Calculates Maxwellian reaction rate for a beam with atomic mass “m_beam”, energy “E_beam”, firing into a target with atomic mass “m_bckg” and temperature “T”.

The “sigma_fun” argument must be a function for a specific charge and n-level of the beam particles. Ref: FIDASIM atomic_tables.f90 bt_maxwellian_n_m.

Parameters

sigma_fun – :py:meth Function to compute a specific cross section [$cm^2$], function of energy/amu ONLY. Expected call form: sigma_fun(erel/ared)
Ti – float, 1D or 2D array Target temperature [keV]. Results will be computed for each Ti value in a vectorized manner.
E_beam – float Beam energy [keV]
m_bckg – float Target atomic mass [amu]
m_beam – float Beam atomic mass [amu]
n_level – int n-level of beam. This is used to evaluate the hydrogen ionization potential, below which an electron is unlikely to charge exchange with surrounding ions.

Returns

output reaction rate in [cm^3/s] units

Return type

rate

aurora.nbi_neutrals.get_NBI_imp_cxr_q(neut_fsa, q, rhop_Ti, times_Ti, Ti_prof, include_fast=True, include_halo=True, debug_plots=False)[source]¶

Compute flux-surface-averaged (FSA) charge exchange recombination for a given impurity with neutral beam components, applying appropriate Maxwellian averaging of cross sections and obtaining rates in [$s^-1$] units. This method expects all neutral components to be given in a dictionary with a structure that is independent of NBI model.

Note that while Ti may be time-dependent, with a time base given by times_Ti, the FSA neutrals are expected to be time-independent. Hence, the resulting CXR rates will only have time dependence that reflects changes in Ti, but not the NBI.

Parameters

neut_fsa – dict Dictionary containing FSA neutral densities in the form that is output by get_neutrals_fsa().
q – int or float Charge of impurity species
rhop_Ti – array-like Sqrt of poloidal flux radial coordinate for Ti profiles.
times_Ti – array-like Time base on which Ti_prof is given [s].
Ti_prof – array-like Ion temperature profile on the rhop_Ti, times_Ti bases.
include_fast – bool, optional If True, include CXR rates from fast NBI neutrals. Default is True.
include_halo – bool, optional If True, include CXR rates from themral NBI halo neutrals. Default is True.
debug_plots – bool, optional If True, plot several plots to assess the quality of the calculation.

Returns

dict

Dictionary containing CXR rates from NBI neutrals. This dictionary has analogous form to the get_neutrals_fsa() function, e.g. we have

rates[beam][f'n={n_level}']['halo']

Rates are on a radial grid corresponding to the input neut_fsa[‘rhop’].

Return type

rates

For details on inputs and outputs, it is recommendeded to look at the internal plotting functions.

aurora.nbi_neutrals.get_ls_cycle()[source]¶

aurora.nbi_neutrals.get_neutrals_fsa(neutrals, geqdsk, debug_plots=True)[source]¶

Compute charge exchange recombination for a given impurity with neutral beam components, obtaining rates in [$s^-1$] units. This method expects all neutral components to be given in a dictionary with a structure that is independent of NBI model (i.e. coming from FIDASIM, NUBEAM, pencil calculations, etc.).

Parameters

neutrals –
dict Dictionary containing fields {“beams”,”names”,”R”,”Z”, beam1, beam2, etc.} Here beam1,beam2,etc. are the names in neutrals[“beams”]. “names” are the names of each beam component, e.g. ‘fdens’,’hdens’,’halo’, etc., ordered according to “names”. “R”,”Z” are the major radius and vertical coordinates [cm] on which neutral density components are given in elements such as
```
neutrals[beams[0]]["n=0"][name_idx]
```
It is currently assumed that n=0,1 and 2 beam components are provided by the user.
geqdsk – gEQDSK post-processed dictionary, as given by the omfit_eqdsk package.
debug_plots – bool, optional If True, various plots are displayed.

Returns

dict

Dictionary of flux-surface-averaged (FSA) neutral densities, in the same units as in the input. Similarly to the input “neutrals”, this dictionary has a structure like

neutrals_ext[beam][f'n={n_level}'][name_idx]

Return type

neut_fsa

aurora.nbi_neutrals.rotation_matrix(alpha, beta, gamma)[source]¶: See the table of all rotation possiblities, on the Tait Bryan side https://en.wikipedia.org/wiki/Euler_angles#Tait.E2.80.93Bryan_angles

aurora.nbi_neutrals.tt_rate_maxwell_average(sigma_fun, Ti, m_i, m_n, n_level)[source]¶

Calculates Maxwellian reaction rate for an interaction between two thermal populations, assumed to be of neutrals (mass m_n) and background ions (mass m_i).

The ‘sigma_fun’ argument must be a function for a specific charge and n-level of the neutral particles. This allows evaluation of atomic rates for charge exchange interactions between thermal beam halos and background ions.

Parameters

sigma_fun – python function Function to compute a specific cross section [cm^2], function of energy/amu ONLY. Expected call form: sigma_fun(erel/ared)
Ti – float or 1D array background ion and halo temperature [keV]
m_i – float mass of background ions [amu]
m_n – float mass of neutrals [amu]
n_level – int n-level of beam. This is used to evaluate the hydrogen ionization potential, below which an electron is unlikely to charge exchange with surrounding ions.
TODO – add effect of toroidal rotation! This will require making the integration in this
2-dimensional. (function) –

Returns

float or 1D array: output reaction rate in [cm^3/s] units

Return type

rate

aurora.nbi_neutrals.uvw_xyz(u, v, w, origin, R)[source]¶

Computes array elements by multiplying the rows of the first array by the columns of the second array. The second array must have the same number of rows as the first array has columns. The resulting array has the same number of rows as the first array and the same number of columns as the second array.

See uvw_to_xyz in fidasim.f90

aurora.nbi_neutrals.xyz_uvw(x, y, z, origin, R)[source]¶

Computes array elements by multiplying the rows of the first array by the columns of the second array. The second array must have the same number of rows as the first array has columns. The resulting array has the same number of rows as the first array and the same number of columns as the second array.

See xyz_to_uvw in fidasim.f90

aurora.janev_smith_rates module¶

Script collecting rates from Janev & Smith, NF 1993. These are useful in aurora to compute total (n-unresolved) charge exchange rates between heavy ions and neutrals.

sciortino, 2020

aurora.janev_smith_rates.js_sigma(E, q, n1, n2=None, type='cx')[source]¶

Cross sections for collisional processes between beam neutrals and highly-charged ions, from Janev & Smith 1993.

Parameters

E – float Normalized beam energy [keV/amu]
q – int Impurity charge before interaction (interacting ion is A^{q+})
n1 – int Principal quantum number of beam hydrogen.
n2 – int Principal quantum number of excited. This may not be needed for some transitions (if so, leave to None).
type – str Type of interaction. Possible choices: {‘exc’,’ioniz’,’cx’} where ‘cx’ refers to electron capture / charge exchange.

Returns

float: Cross section of selected process, in [cm^2] units.

Return type

sigma

See comments in Janev & Smith 1993 for uncertainty estimates.

aurora.janev_smith_rates.js_sigma_cx_n1_q1(E)[source]¶: Electron capture cross section for H^{+} + H(1s) –> H + H^+ Section 2.3.1

aurora.janev_smith_rates.js_sigma_cx_n1_q2(E)[source]¶: Electron capture cross section for He^{2+} + H(1s) –> He^+ + H^+ Section 3.3.1

aurora.janev_smith_rates.js_sigma_cx_n1_q4(E)[source]¶: Electron capture cross section for Be^{4+} + H(1s) –> Be^{3+} + H^+ Section 4.3.1

aurora.janev_smith_rates.js_sigma_cx_n1_q5(E)[source]¶: Electron capture cross section for B^{5+} + H(1s) –> B^{4+} + H^+ Section 4.3.2

aurora.janev_smith_rates.js_sigma_cx_n1_q6(E)[source]¶: Electron capture cross section for C^{6+} + H(1s) –> C^{5+} + H^+ Section 4.3.3

aurora.janev_smith_rates.js_sigma_cx_n1_q8(E)[source]¶: Electron capture cross section for O^{8+} + H(1s) –> O^{7+} + H^+ Section 4.3.4

aurora.janev_smith_rates.js_sigma_cx_n1_qg8(E, q)[source]¶: Electron capture cross section for A^{q+} + H(1s) –> A^{(q-1)+} + H^+, q>8 Section 4.3.5, p.172

aurora.janev_smith_rates.js_sigma_cx_n2_q2(E)[source]¶: Electron capture cross section for He^{2+} + H(n=2) –> He^+ + H^+ Section 3.3.2

aurora.janev_smith_rates.js_sigma_cx_ng1_q1(E, n1)[source]¶: Electron capture cross section for H^{+} + H(n) –> H + H^+ , n>1 Section 2.3.2

aurora.janev_smith_rates.js_sigma_cx_ng1_qg3(E, n1, q)[source]¶: Electron capture cross section for A^{q+} + H^*(n) –> A^{(q-1)+}+H^+ , q>3 Section 4.3.6, p.174

aurora.janev_smith_rates.js_sigma_cx_ng2_q2(E, n1)[source]¶: Electron capture cross section for He^{2+} + H*(n) –> He^+ + H^+ , n>2 Section 3.2.3

aurora.janev_smith_rates.js_sigma_ioniz_n1_q8(E)[source]¶: Ionization cross section for O^{8+} + H(1s) –> O^{8+} + H^+ +e^- Section 4.2.4

aurora.janev_smith_rates.plot_js_sigma(q=18)[source]¶

aurora.synth_diags module¶

aurora.synth_diags.centrifugal_asym(rhop, Rlfs, omega, Zeff, A_imp, Z_imp, Te, Ti, main_ion_A=2, plot=False, nz=None, geqdsk=None)[source]¶

Estimate impurity poloidal asymmetry effects from centrifugal forces.

The result of this function is $\lambda$, defined such that

\[n(r, heta) = n_0(r) imes \exp\left{\lambda(\]

ho) (R(r, heta)^2- R_0^2) ight}

See Odstrcil et al. 2018 Plasma Phys. Control. Fusion 60 014003 for details on centrifugal asymmetries. Also see Appendix A of Angioni et al 2014 Nucl. Fusion 54 083028 for details on these should also be accounted for when comparing transport coefficients used in Aurora (on a rvol grid) to coefficients used in codes that use other coordinate systems (e.g. based on rmid).

Args:

rhoparray (nr,)
Sqrt of normalized poloidal flux grid.

Rlfsarray (nr,)
Major radius on the Low Field Side (LFS), at points corresponding to rhop values

omegaarray (nt,nr) or (nr,) [ rad/s ]
Toroidal rotation on Aurora temporal time_grid and radial rhop_grid (or, equivalently, rvol_grid) grids.

Zeffarray (nt,nr), (nr,) or float
Effective plasma charge on Aurora temporal time_grid and radial rhop_grid (or, equivalently, rvol_grid) grids. Alternatively, users may give Zeff as a float (taken constant over time and space).

A_impfloat
Impurity ion atomic mass number (e.g. 40 for Ca)

Z_imparray (nr, ) or int

Charge state of the impurity of interest. This can be an array, giving the expected charge state at every radial position, or just a float.

Tearray (nr,nt)
Electron temperature (eV)

Tiarray (nr, nt)
Background ion temperature (eV)

main_ion_Aint, optional
Background ion atomic mass number. Default is 2 for D.

Keyword Args:

plotbool
If True, plot asymmetry factor $\lambda$ vs. radius and show the predicted 2D impurity density distribution at the last time point.

nzarray (nr,nZ)
Impurity charge state densities (output of Aurora at a specific time slice), only used for 2D plotting.

geqdskdict
Dictionary containing the omfit_eqdsk reading of the EFIT g-file.

Returns:

CF_lamarray (nr,)
Asymmetry factor, defined as $\lambda$ in the expression above.

aurora.synth_diags.line_int_weights(R_path, Z_path, rhop_path, dist_path, R_axis=None, rhop_out=None, CF_lam=None)[source]¶

Obtain weights for line integration on a rhop grid, given the 3D path of line integration in the (R,Z,Phi) coordinates, as well as the value of sqrt of normalized poloidal flux at each point along the path.

Parameters

R_path – array (np,) Values of the R coordinate [m] along the line integration path.
Z_path – array (np,) Values of the Z coordinate [m] along the line integration path.
rhop_path – array (np,) Values of the rhop coordinate along the line integration path.
dist_path – array (np,) Vector starting from 0 to the maximum distance [m] considered along the line integration.

Keyword Arguments

R_axis – float R value at the magnetic axis [m]. Only used for centrifugal asymmetry effects if CF_lam is not None.
rhop_out – array (nr,) The sqrt of normalized poloidal flux grid on which weights should be computed. If left to None, an equally-spaced grid with 201 points from the magnetic axis to the LCFS is used.
CF_lam – array (nr,) Centrifugal (CF) asymmetry exponential factor, returned by the centrifugal_asym() function. If provided, this is taken to be on an input rhop_out grid. If left to None, no CF asymmetry is considered.

Module contents¶

aurora