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)

Bases: object

Setup the input dictionary for an Aurora ion transport simulation from the given namelist.

Parameters:

• namelist (dict) – Dictionary containing aurora inputs. See default_nml.py for some defaults, which users should modify for their runs.

• geqdsk (dict, optional) – EFIT gfile as returned after postprocessing by the omfit_classes.omfit_eqdsk package (OMFITgeqdsk class). If left to None (default), the minor and major radius must be indicated in the namelist in order to create a radial grid.

calc_Zeff()

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)

Estimate impurity poloidal asymmetry effects from centrifugal forces. See notes the centrifugal_asymmetry() 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).

• plot (bool) – If True, plot asymmetry factor \(\lambda\) vs. radius

Returns:

CF_lambda – Asymmetry factor, defined as \(\lambda\) in the centrifugal_asymmetry() function docstring.

Return type:

array (nr,)

check_conservation(plot=True, axs=None, plot_resolutions=False)

Check particle conservation for an aurora simulation.

Parameters:

• plot (bool, optional) – If True, plot time histories in each particle reservoir and display quality of particle conservation.

• 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:

• out (dict) – Dictionary containing density of particles in each reservoir.

• axs (matplotlib.Axes instances, only returned if plot=True) – Array-like structure containing two matplotlib.Axes instances, (ax1,ax2). See optional input argument.

get_aurora_kin_profs(min_T=1.01, min_ne=10000000000.0)

Get kinetic profiles on radial and time grids.

get_par_loss_rate(trust_SOL_Ti=False)

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

Parameters:

trust_SOL_Ti (bool) – If True, the input Ti is trusted also in the SOL to calculate a parallel loss rate. Often, Ti measurements in the SOL are unrealiable, so this parameter is set to False by default.

Returns:

dv – Parallel loss rates in \(s^{-1}\) units. Values are zero in the core region and non-zero in the SOL.

Return type:

array (space,time)

interp_kin_prof(prof)

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.

load(filename)

Load aurora_sim object.

load_dict(aurora_dict)

plot_resolutions()

Convenience function to show time and spatial resolution in Aurora simulation setup.

reload_namelist(namelist=None)

(Re-)load namelist to update scalar variables.

run_aurora(D_z, V_z, times_DV=None, nz_init=None, unstage=True, alg_opt=1, evolneut=False, use_julia=False, plot=False)

Run a simulation using the provided diffusion and convection profiles as a function of space, time and potentially also ionization state. Users can 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 (array, shape of (space,time,nZ) or (space,time) or (space,)) – Diffusion coefficients, in units of \(cm^2/s\). This may be given as a function of space only, (space,time) or (space,nZ, time), where nZ indicates the number of charge states. If given with 1 or 2 dimensions, it is assumed that all charge states should have the same diffusion coefficients. If given as 1D, it is further assumed that diffusion is time-independent. Note that it is assumed that radial profiles are already on the self.rvol_grid radial grid.

• V_z (array, shape of (space,time,nZ) or (space,time) or (space,)) – Convection coefficients, in units of \(cm/s\). This may be given as a function of space only, (space,time) or (space,nZ, time), where nZ indicates the number of charge states. If given with 1 or 2 dimensions, it is assumed that all charge states should have the same convection coefficients. If given as 1D, it is further assumed that convection is time-independent. Note that it is assumed that radial profiles are already on the self.rvol_grid radial grid.

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

• unstage (bool, optional) – If superstages are indicated in the namelist, this parameter sets whether the output should be “unstaged” by multiplying by the appropriate fractional abundances of all charge states at ionization equilibrium. Note that this unstaging process cannot account for transport and is therefore only an approximation, to be used carefully.

• alg_opt (int, optional) – If alg_opt=1, use the finite-volume algorithm proposed by Linder et al. NF 2020. If alg_opt=0, 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. NB: It is recommended to only use this with explicit 2D sources, otherwise

• 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:

• nz (array, (nr,nZ,nt)) – Charge state densities [\(cm^{-3}\)] over the space and time grids. If a number of superstages are indicated in the input, only charge state densities for these are returned.

• N_wall (array (nt,)) – Number of particles at the wall reservoir over time.

• N_div (array (nt,)) – Number of particles in the divertor reservoir over time.

• N_pump (array (nt,)) – Number of particles in the pump reservoir over time.

• N_ret (array (nt,)) – Number of particles temporarily held in the wall reservoirs.

• N_tsu (array (nt,)) – Edge particle loss [\(cm^{-3}\)]

• N_dsu (array (nt,)) – Parallel particle loss [\(cm^{-3}\)]

• N_dsul (array (nt,)) – Parallel particle loss at the limiter [\(cm^{-3}\)]

• rcld_rate (array (nt,)) – Recycling from the divertor [\(s^{-1} cm^{-3}\)]

• rclw_rate (array (nt,)) – Recycling from the wall [\(s^{-1} cm^{-3}\)]

run_aurora_steady(D_z, V_z, nz_init=None, unstage=False, alg_opt=1, evolneut=False, use_julia=False, tolerance=0.01, max_sim_time=100, dt=0.0001, dt_increase=1.05, n_steps=100, plot=False)

Run an Aurora simulation until reaching steady state profiles. This method calls run_aurora() checking at every iteration whether profile shapes are still changing within a given fractional tolerance. Note that this method differs from run_aurora_steady_analytic() in that it runs time-dependent simulations until reaching steady-state, rather than using an analytic solution for steady-state profiles.

Parameters:

• D_z (array, shape of (space,nZ) or (space,)) – Diffusion coefficients, in units of \(cm^2/s\). This may be given as a function of space only or (space,nZ). No time dependence is allowed in this function. Here, nZ indicates the number of charge states. Note that it is assumed that radial profiles are already on the self.rvol_grid radial grid.

• V_z (array, shape of (space,nZ) or (space,)) – Convection coefficients, in units of \(cm/s\). This may be given as a function of space only or (space,nZ). No time dependence is allowed in this function. Here, nZ indicates the number of charge states.

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

• unstage (bool, optional) – If a list of superstages are provided in the namelist, this parameter sets whether the output should be “unstaged”. See docs for run_aurora() for details.

• alg_opt (int, optional) – If alg_opt=1, use the finite-volume algorithm proposed by Linder et al. NF 2020. If alg_opt=0, 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. See docs for run_aurora() for details.

• use_julia (bool, optional) – If True, run the Julia pre-compiled version of the code. See docs for run_aurora() for details.

• tolerance (float) – Fractional tolerance in charge state profile shapes. This method reports charge state density profiles obtained when the discrepancy between normalized profiles at adjacent time steps varies by less than this tolerance fraction.

• max_sim_time (float) – Maximum time in units of seconds for which simulations should be run if a steady state is not found.

• dt (float) – Initial time step to apply, in units of seconds. This can be increased by a multiplier given by :param:`dt_increase` after each time step.

• dt_increase (float) – Multiplier for time steps.

• n_steps (int) – Number of time steps (>2) before convergence is checked.

• plot (bool) – If True, plot time evolution of charge state density profiles to show convergence.

run_aurora_steady_analytic(D_z, V_z)

Evaluate the analytic steady state solution of the transport equation. Small differences in absolute densities from the full time-dependent, multi-reservoir Aurora solutions can be caused by recycling and divertor models, which are not included here.

Parameters:

• D_z (array, shape of (space,nZ) or (space,)) – Diffusion coefficients, in units of \(cm^2/s\). This may be given as a function of space only or (space,nZ). No time dependence is allowed in this function. Here, nZ indicates the number of charge states. Note that it is assumed that radial profiles are already on the self.rvol_grid radial grid.

• V_z (array, shape of (space,nZ) or (space,)) – Convection coefficients, in units of \(cm/s\). This may be given as a function of space only or (space,nZ). No time dependence is allowed in this function. Here, nZ indicates the number of charge states.

save(filename)

Save state of aurora_sim object.

save_dict()

set_time_dept_atomic_rates(superstages=[], metastables=False)

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.

Parameters:

• superstages (list or 1D array) – Indices of charge states that should be kept as superstages. The default is to have this as an empty list, in which case all charge states are kept.

• metastables (bool) – Load metastable resolved atomic data and rates. Default is False.

Sne_rates

Effective ionization rates [s]. If superstages were indicated, these are the rates of superstages.

Type:

array (space, nZ(-super), time)

Rne_rates

Effective recombination rates [s]. If superstages were indicated, these are the rates of superstages.

Type:

array (space, nZ(-super), time)

Qne_rates

Cross-coupling coefficients, only computed if :param:`metastables`=True.

Type:

array (space, nZ(-super), time)

Xne_rates

Parent cross-coupling coefficients, only computed if :param:`metastables`=True.

Type:

array (space, nZ(-super), time)

setup_grids()

Method to set up radial and temporal grids given namelist inputs.

setup_kin_profs_depts()

Method to set up Aurora inputs related to the kinetic background from namelist inputs.

superstage_DV(D_z, V_z, times_DV=None, opt=1)

Reduce the dimensionality of D and V time-dependent profiles for the case in which superstaging is applied.

Three options are currently available:

1. opt=1 gives a simple selection of D_z and V_z fields corresponding to each superstage index.

2. opt=2 averages D_z and V_z over the charge states that are part of each superstage.

#. opt=3 weights D_z and V_z corresponding to each superstage by the fractional abundances at ionization equilibrium. This is mostly untested – use with care!

Parameters:

• D_z (array, shape of (space,time,nZ)) – Diffusion coefficients, in units of \(cm^2/s\).

• V_z (array, shape of (space,time,nZ)) – Convection coefficients, in units of \(cm/s\).

Returns:

• Dzf (array, shape of (space,time,nZ-superstages)) – Diffusion coefficients of superstages, in units of \(cm^2/s\).

• Vzf (array, shape of (space,time,nZ-superstages)) – Convection coefficients of superstages, in units of \(cm/s\).

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)

Bases: object

Fast linear interpolation for 1D and 2D vector data on equally spaced grids. This offers optimal speed in Python for interpolation of atomic data tables such as the ADAS ones.

Parameters:

• grids (list of arrays, N=len(grids), N=1 or N=2) – List of 1D arrays with equally spaced grid values for each dimension

• values (N+1 dimensional array of values used for interpolation) – Values to interpolate. The first dimension typically refers to different ion stages, for which data is provided on the input grids. Other dimensions refer to values on the density and temperature grids.

class aurora.atomic.adas_file(filepath)

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.

Parameters:

filepath (str) – Path to location where ADAS file is located.

load()

ADF11 format description: https://www.adas.ac.uk/man/appxa-11.pdf

plot(fig=None, axes=None)

Plot data from input ADAS file. If provided, the arguments allow users to overplot and compare data from multiple files.

Parameters:

• fig (matplotlib Figure object) – If provided, add specification as to which ADAS file is being plotted.

• axes (matplotlib Axes object (or equivalent)) – If provided, plot on these axes. Note that this typically needs to be a set of axes for each plotted charge state. Users may want to call this function once first to get some axes, and then pass those same axes to a second call for another file to compare with.

aurora.atomic.get_adas_file_types()

Obtain a description of each ADAS file type and its meaning in the context of Aurora.

Returns:

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

Return type:

dict

Notes

For background on ADAS generalized collisional-radiative modeling and data formats, refer to [1]_.

References

[1]

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

aurora.atomic.get_atom_data(imp, files=['acd', 'scd'])

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.

• files (list or dict) –

  ADAS file types to be fetched. Default is [“acd”,”scd”] for effective ionization and recombination rates (excluding CX) using default files, listed in adas_files_adas_files_dict(). If users prefer to use specific files, they may pass a dictionary instead, of the form

  {'acd': 'acd89_ar.dat', 'scd': 'scd89_ar.dat'}
  

  or

  {'acd': 'acd89_ar.dat', 'scd': None}
  

  if only some of the files need specifications and others (given as None) should be taken from the default files.

Returns:

atom_data – 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 res.logNe, res.logT, res.logdata, res.meta_ind, res.metastables

Return type:

dict

aurora.atomic.get_atomic_relax_time(atom_data, ne_cm3, Te_eV=None, Ti_eV=None, n0_by_ne=0.0, superstages=[], tau_s=inf, plot=True, ax=None, ls='-')

Obtain the relaxation time of the ionization equilibrium for a given atomic species.

If n0_by_ne is not 0, thermal charge exchange is added to radiative and dielectronic recombination. This function can work with ne,Te and n0_by_ne arrays of arbitrary dimension. It uses a matrix SVD approach in order to find the relaxation rates, as opposed to the simpler approach of get_frac_abundances(), but fractional abundances produced by the two methods should always be the same.

This function also allows use of superstages as well as specification of a \(\tau\) value representing the effect of transport on the ionization equilibrium. NB: this is only a rough metric to characterize complex physics.

Parameters:

• atom_data (dictionary of atomic ADAS files (only acd, scd are required; ccd is) – necessary only if n0_by_ne is not 0).

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

• Ti_eV (float or array) – Bulk ion temperature in units of eV, only needed for CX. If left to None, Ti is set equal to Te.

• n0_by_ne (float or array, optional) – Ratio of background neutral hydrogen to electron density. If set to 0, CX is not considered.

• superstages (list or 1D array) – Indices of charge states of chosen ion that should be included. If left empty, all ion stages are included. If only some indices are given, these are modeled as “superstages”.

• tau_s (float, opt) – Value of the particle residence time [s]. This is a scalar value that can be used to model the effect of transport on ionization equilibrium. Setting tau=np.inf (default) corresponds to no effect from transport.

• plot (bool, optional) – If True, the atomic relaxation time is plotted as a function of Te. Default is True.

• ax (matplotlib.pyplot Axes instance) – Axes on which to plot if plot=True. If False, new axes are created.

• ls (str, optional) – Line style for plots. Continuous lines are used by default.

Returns:

• Te (array) – electron temperatures as a function of which the fractional abundances and rate coefficients are given.

• fz (array, (space,nZ)) – Fractional abundances across the same grid used by the input ne,Te values.

• rate_coeffs (array, (space, nZ)) – Rate coefficients in units of [\(s^{-1}\)].

Examples

To visualize relaxation times for a given species: >>> atom_data = aurora.atomic.get_atom_data(‘N’, [“scd”, “acd”]) >>> aurora.get_atomic_relax_time(atom_data, [1e14], Te_eV=np.linspace(0.1,200,1000), plot=True)

To compare ionization balance with different values of ne*tau: >>> Te0, fz0, r0 = aurora.get_atomic_relax_time(atom_data, [1e14], Te_eV=np.linspace(0.1,200,1000), tau_s=1e-3, plot=False) >>> Te1, fz1, r1 = aurora.get_atomic_relax_time(atom_data, [1e14], Te_eV=np.linspace(0.1,200,1000), tau_s=1e-2, plot=False) >>> Te2, fz2, r2 = aurora.get_atomic_relax_time(atom_data, [1e14], Te_eV=np.linspace(0.1,200,1000), tau_s=1e-1, plot=False) >>> >>> plt.figure() >>> for cs in np.arange(fz0.shape[1]): >>> l = plt.plot(Te0, fz0[:,cs], ls=’-‘) >>> plt.plot(Te1, fz1[:,cs], c=l[0].get_color(), ls=’–‘) >>> plt.plot(Te2, fz2[:,cs], c=l[0].get_color(), ls=’-.’)

aurora.atomic.get_cs_balance_terms(atom_data, ne_cm3=50000000000000.0, Te_eV=None, Ti_eV=None, include_cx=True, metastables=False)

Get S*ne, R*ne and cx*ne rates 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.

• Ti_eV (float or array) – Bulk ion temperature in units of eV, only needed for CX. If left to None, Ti is set equal to Te.

• include_cx (bool) – If True, obtain charge exchange terms as well.

Returns:

• Te (array (n_Te)) – Te grid on which atomic rates are given

• Sne, Rne (,cxne, Qne, Xne) (arrays (n_ne,n_Te)) – atomic rates for effective ionization, radiative+dielectronic recombination (+ charge exchange, + crosscoupling if requested). All terms will be in units of \(s^{-1}\).

Notes

The nomenclature with ‘ne’ at the end of rate names indicates that these are rates in units of \(m^3\cdot s^{-1}\) multiplied by the electron density ‘ne’ (hence, final units of \(s^{-1}\)).

aurora.atomic.get_frac_abundances(atom_data, ne_cm3, Te_eV=None, Ti_eV=None, n0_by_ne=0.0, superstages=[], plot=True, ax=None, rho=None, rho_lbl=None)

Calculate fractional abundances from ionization and recombination equilibrium. If n0_by_ne is not 0, 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.

• Ti_eV (float or array, optional) – Bulk ion temperature in units of eV. If left to None, Ti is set to be equal to Te

• n0_by_ne (float or array, optional) – Ratio of background neutral hydrogen to electron density. If not 0, CX is considered.

• superstages (list or 1D array) – Indices of charge states of chosen ion that should be included. If left empty, all ion stages are included. If only some indices are given, these are modeled as “superstages”.

• 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 “x”.

Returns:

• Te (array) – electron temperatures as a function of which the fractional abundances and rate coefficients are given.

• fz (array, (space,nZ)) – Fractional abundances across the same grid used by the input ne,Te values.

aurora.atomic.get_natural_partition(ion, plot=True)

Identify natural partition of charge states by plotting the variation of ionization energy as a function of charge for a given ion, using the ADAS metric \(2 (I_{z+1}-I_z)/(I_{z+1}+I_z)\).

Parameters:

• ion (str) – Atomic symbol of species of interest.

• plot (bool) – If True, plot the variation of ionization energy.

Returns:

q – Metric to identify natural partition.

Return type:

1D array

Notes

A ColRadPy installation must be available for this function to work.

aurora.atomic.gff_mean(Z, Te)

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, freq='all', cutoff=0.1)

Approximate impurity bremsstrahlung in :math:`W/m^3`for a given range of frequencies or a specific frequency.

We apply here a cutoff for Bremsstrahlung at h*c/lambda = cutoff*Te, where cutoff is an input parameter, conventionally set to 0.1 (default).

Gaunt factors from gff_mean() are applied. NB: recombination is not included.

Formulation based on Hutchinson’s Principles of Plasma Diagnostics, p. 196, Eq. (5.3.40).

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}\)]

• freq (float, 1D array, or str) – If a float, calculate bremsstrahlung from all charge states at this frequency. If a 1D array, evaluate bremsstrahlung at these wavelengths. If set to all, then bremsstrahlung is integrated over the whole range from plasma frequency to cutoff Frequencies are expected in units of \(s^{-1}\).

• cutoff (float) – Fraction of Te below which bremsstrahlung is set to 0. A value of 0.1 is commonly set and is the default.

Returns:

Bremsstrahlung for each charge state at the given frequency or, if multiple frequences are given (or if freq=’all’), integrated over frequencies. Units of \(W/cm^3\).

Return type:

array (time,nZ,space)

aurora.atomic.interp_atom_prof(atom_table, xprof, yprof, log_val=False, x_multiply=True)

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

Parameters:

• atom_table (list) – object atom_data, containing atomic data from one of the ADAS files.

• xprof (array (nt,nr)) – Spatio-temporal profiles of the first coordinate of the ADAS file table (usually electron density in \(cm^{-3}\))

• yprof (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^{xprof}\).

Returns:

interp_vals – Interpolated atomic data on time,charge state and spatial grid that correspond to the ion of interest and the spatiotemporal grids of xprof and yprof.

Return type:

array (nt,nion,nr)

Notes

This function uses np.log10 and exponential operations to optimize speed, since it has been observed that base-e operations are faster than base-10 operations in numpy.

aurora.atomic.null_space(A)

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)

Compare effective ionization rate for each charge state with the characteristic transit time that passing and trapped impurity ions take to travel a parallel distance \(L = q R\), defining

\[\nu_{ion}^* \equiv \nu_{ion} \tau_t = \nu_{ion} \frac{q R}{v_{th}} = \frac{\sum_z n_z \nu_z^{ion}}{\sum_z n_z} q R \sqrt{\frac{m_{imp}}{2 k_B T_i}}\]

following Ref.[1]_. If the normalized ionization rate (\(\nu_{ion}^*\)) is less than 1, then flux surface averaging of background asymmetries (e.g. from edge or beam neutrals) may be taken as a good approximation of reality; in this case, 1.5D simulations of impurity transport are expected to be valid. If, on the other hand, \(\nu_{ion}^*>1\) then local effects may be too important to ignore.

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:

nu_ioniz_star – 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:

array (r,cs) or (r,)

References

[1]

R.Dux et al. Nucl. Fusion 60 (2020) 126039

aurora.atomic.read_adf12(filename, block, Ebeam, ne_cm3, Ti_eV, zeff)

Read charge exchange effective emission coefficients from ADAS ADF12 files.

Files may be automatically downloaded using :py:fun:`~aurora.adas_files.get_adas_file_loc`.

Parameters:

• filename (str) – adf12 file name/path

• block (int) – Source block selected

• Ebeam (float) – Energy of the neutral beam population of interest, in units of \(eV/amu\).

• ne_cm3 (float or 1D array) – Electron densities at which to evaluate coefficients, in units of \(cm^{-3}\).

• Ti_eV (float or 1D array) – Bulk ion temperature at which to evaluate coefficients, in units of \(eV\).

• Zeff (float or 1D array) – Effective background charge.

Returns:

Interpolated coefficients, units of \(cm^3/s\).

Return type:

float or 1D array

aurora.atomic.read_adf21(filename, Ebeam, ne_cm3, Te_eV)

Read ADAS ADF21 or ADF22 files.

ADF21 files contain effective beam stopping/excitation coefficients. ADF22 contain effective beam emission/population coefficients.

Files may be automatically downloaded using :py:fun:`~aurora.adas_files.get_adas_file_loc`.

Parameters:

• filename (str) – adf21 or adf22 file name/path

• Ebeam (float) – Energy of the neutral beam, in units of \(eV/amu\).

• ne_cm3 (float or 1D array) – Electron densities at which to evaluate coefficients.

• Te_eV (float or 1D array) – Electron temperature at which to evaluate coefficients.

Returns:

Interpolated coefficients. For ADF21 files, these have units of :math:`cm^3/s`for ADF21 files. For ADF22, they correspond to n=2 fractional abundances.

Return type:

float or 1D array

aurora.atomic.read_filter_response(filepath, adas_format=True, plot=False)

Read a filter response function over energy.

This function attempts to read the data checking for the following formats (in this order):

1. The ADAS format. Typically, this data is from obtained from http://xray.uu.se and produced via ADAS routines.

2. The format returned by the Center for X-Ray Optics website .

Note that filter response functions are typically a combination of a filter transmissivity and a detector absorption.

Parameters:

• filepath (str) – Path to filter file of interest.

• plot (bool) – If True, the filter response function is plotted.

aurora.atomic.superstage_rates(R, S, superstages, save_time=None)

Compute rate for a set of ion superstages. Input and output rates are log-values in arbitrary base.

Parameters:

• R (array (time,nZ,space)) – Array containing the effective recombination rates for all ion stages, These are typically combinations of radiative and dielectronic recombination, possibly also of charge exchange recombination.

• S (array (time,nZ,space)) – Array containing the effective ionization rates for all ion stages.

• superstages (list or 1D array) – Indices of charge states of chosen ion that should be included.

• save_time (list or 1D array of bools) – Indices of the timeslcies which are actually returned by AURORA

Returns:

• superstages (array) – Set of superstages including 0,1 and final stages if these were missing in the input.

• R_rates_super (array (time,nZ-super,space)) – effective recombination rates for superstages

• S_rates_super (array (time,nZ-super,space)) – effective ionization rates for superstages

• fz_upstage (array (space, nZ, time_) – fractional abundances of stages within superstages

aurora.adas_files module

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

aurora.adas_files.adas_files_dict()

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:

files – 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:

dict

aurora.adas_files.get_adas_file_loc(filename, filetype='adf11')

Find location of requested atomic data file for the indicated ion. Accepts all ADAS “derived” data file types (adf11, adf12, adf13, adf15, adf21, adf22). The search proceeds with the following attempts, in this order:

1. If the file is available in Aurora/adas_data/filetype, with filetype given by the user, always use this data.

2. If the input filename is actually a full path to the file on the local system, use this file and copy it to Aurora/aurora/adas_data/actual_filename, where actual_filename is the file name rather than the full path.

3. If the environmental variable “AURORA_ADAS_DIR” is defined, attempt to find the file there, under adf11/ and adf15/ directories. AURORA_ADAS_DIR may for example be the path to a central repository of ADAS files for Aurora on a cluster where not everyone may have write-permissions. For this option, files are not copied at all.

4. Attempt to fetch the file remotely via open.adas.ac.uk and save it in Aurora/aurora/adas_data/filetype/.

Parameters:

• filename (str) – Name of the ADAS file of interest, e.g. ‘plt89_ar.dat’.

• filetype (str) – ADAS file type. Only derived data types are allowed, i.e. one of “adf11”, “adf12”, “adf13”, “adf15”, “adf21”, “adf22”

Returns:

file_loc – Full path to the requested file.

Return type:

str

aurora.radiation module

aurora.radiation.adf04_files()

Collection of trust-worthy ADAS ADF04 files. This function will be moved and expanded in ColRadPy in the near future.

aurora.radiation.adf15_line_identification(pec_files, lines=None, Te_eV=1000.0, ne_cm3=50000000000000.0, mult=[])

Display all photon emissivity coefficients from the given list of ADF15 files and (optionally) compare to a set of chosen wavelengths, given in units of Angstrom.

Parameters:

• pec_files (str or list of str) – Path to a single ADF15 file or a list of files.

• lines (dict, list or 1D array) – Lines to overplot with the loaded PECs to consider overlap within spectrum. This argument may be a dictionary, with keys corresponding to line names and values corresponding to wavelengths (in units of Angstrom). If provided as a list or array, this is assumed to contain only wavelengths in Angstrom.

• Te_eV (float) – Single value of electron temperature at which PECs should be evaluated [\(eV\)].

• ne_cm3 (float) – Single value of electron density at which PECs should be evaluated [\(cm^{-3}\)].

• mult (list or array) – Multiplier to apply to lines from each PEC file. This could be used for example to rescale the results of multiple ADF15 files by the expected fractional abundance or density of each element/charge state.

Notes

To attempt identification of spectral lines, one can load a set of ADF15 files, calculate approximate fractional abundances at equilibrium and overplot expected emissivities in a few steps:

>>> pec_files = ['mypecs1','mypecs2','mypecs3']
>>> Te_eV=500; ne_cm3=5e13; ion='Ar'   # examples
>>> atom_data = aurora.atomic.get_atom_data(ion,['scd','acd'])
>>> _Te, fz = aurora.atomic.get_frac_abundances(atom_data, ne_cm3, Te_eV, plot=False)
>>> mult = [fz[0,10], fz[0,11], fz[0,12]] # to select charge states 11+, 12+ and 13+, for example
>>> aurora.adf15_line_identification(pec_files, Te_eV=Te_eV, ne_cm3=ne_cm3, mult=mult)

Minimal Working Example

>>> Te_eV=500; ne_cm3=5e13 # eV and cm^{-3}
>>> filepath = aurora.get_adas_file_loc('pec96#he_pju#he0.dat', filetype='adf15')
>>> aurora.adf15_line_identification(filepath, Te_eV=Te_eV, ne_cm3=ne_cm3)

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

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.

• 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:

res – Dictionary containing radiation terms, depending on the activated flags.

Return type:

dict

Notes

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.

aurora.radiation.get_colradpy_pec_prof(ion, cs, rhop, ne_cm3, Te_eV, lam_nm, lam_width_nm, adf04_loc, meta_idxs=[0], pec_threshold=1e-20, pec_units=2, plot=True)

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 ‘17’, indicating total charge of ion (e.g. ‘17’ would be for Li-like Ca)

• rhop (array (nr,)) – Sqrt 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]

• adf04_loc (str) – Location from which ADF04 files listed in adf04_files() should be fetched.

• meta_idxs (list of integers) – List of levels in ADF04 file to be treated as metastable states. Default is [0] (only ground state).

• prec_threshold (float) – Minimum value of PECs to be considered, in \(photons \cdot cm^3/s\)

• pec_units (int) – If 1, results are given in \(photons \cdot cm^3/s\); if 2, they are given in \(W.cm^3\). Default is 2.

• plot (bool) – If True, plot lines profiles and total.

Returns:

pec_tot_prof – Radial profile of PEC intensity, in units of \(photons \cdot cm^3/s\) (if pec_units=1) or \(W \cdot cm^3\) (if pec_units=2).

Return type:

array (nr,)

aurora.radiation.get_cooling_factors(imp, ne_cm3, Te_eV, n0_cm3=0.0, ion_resolved=False, superstages=[], line_rad_file=None, cont_rad_file=None, sxr=False, plot=True, ax=None)

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

Parameters:

• imp (str) – Atomic symbol of ion of interest

• ne_cm3 (1D array) – Electron density [\(cm^{-3}\)], used to find charge state fractions at ionization equilibrium.

• Te_eV (1D array) – Electron temperature [\(eV\)] at which cooling factors should be obtained.

• n0_cm3 (1D array or float) – Background H/D/T neutral density [\(cm^{-3}\)] used to account for charge exchange when calculating ionization equilibrium. If left to 0, charge exchange effects are not included.

• ion_resolved (bool) – If True, cooling factors are returned for each charge state. If False, they are summed over charge states. The latter option is useful for modeling where charge states are assumed to be in ionization equilibrium (no transport). Default is False.

• superstages (list or 1D array) – List of superstages to consider. An empty list (default) corresponds to the inclusion of all charge states. Note that when ion_resolved=False, cooling coefficients are independent of whether superstages are being used or not.

• line_rad_file (str or None) – Location of ADAS ADF11 file containing line radiation data. This can be a PLT (unfiltered) or PLS (filtered) file. If left to None, the default file given in adas_files_dict() will be used.

• cont_rad_file (str or None) – Location of ADAS ADF11 file containing recombination and bremsstrahlung radiation data. This can be a PRB (unfiltered) or PRS (filtered) file. If left to None, the default file given in adas_files_dict() will be used.

• sxr (bool) – If True, line radiation, recombination and bremsstrahlung radiation are taken to be from SXR-filtered ADAS ADF11 files, rather than from unfiltered files.

• plot (bool) – If True, plot all radiation components, summed over charge states.

• ax (matplotlib.Axes instance) – If provided, plot results on these axes.

Returns:

• line_rad_tot (1D array) – Cooling coefficient from line radiation [\(W\cdot m^3\)]. Depending on whether :param:`sxr`=True or False, this indicates filtered or unfiltered radiation, respectively.

• cont_rad_tot (1D array) – Cooling coefficient from continuum radiation [\(W\cdot m^3\)]. Depending on whether `sxr`=True or False, this indicates filtered or unfiltered radiation, respectively.

aurora.radiation.get_local_spectrum(adf15_file, ne_cm3, Te_eV, ion_exc_rec_dens=[0, 1, 0], Ti_eV=None, n0_cm3=0.0, dlam_A=0.0, plot_spec_tot=True, plot_all_lines=False, no_leg=False, ax=None)

Plot spectrum based on the lines contained in an ADAS ADF15 file at specific values of electron density and temperature. Charge state densities can be given explicitely, or alternatively charge state fractions will be automatically computed from ionization equilibrium (no transport).

Parameters:

• adf15_file (str or dict) – Path on disk to the ADAS ADF15 file of interest or dictionary returned when calling the read_adf15() with this file path as an argument. All wavelengths and radiating components in the file or dictionary will be read/processed.

• ne_cm3 (float) – Local value of electron density, in units of \(cm^{-3}\).

• Te_eV (float) – Local value of electron temperature, in units of \(eV\). This is used to evaluate local values of photon emissivity coefficients.

• ion_exc_rec_dens (list of 3 floats) – Density of ionizing, excited and recombining charge states that may contribute to emission from the given ADF15 file. In the absence of charge state densities from particle transport modeling, these scalars may be taken from the output of aurora.atomic.get_frac_abundances(). Default is [0,1,0], which means that only excitation components are considered.

• Ti_eV (float) – Local value of ion temperature, in units of \(eV\). This is used to represent the effect of Doppler broadening. If left to None, it is internally set equal to Te_eV.

• n0_cm3 (float, optional) – Local density of atomic neutral hydrogen isotopes. This is only used if the provided ADF15 file contains charge exchange contributions.

• dlam_A (float or 1D array) – Doppler shift in A. This can either be a scalar or an array of the same shape as the output wavelength array. For the latter option, it is recommended to call this function twice to find the wave_final_A array first.

• plot_spec_tot (bool) – If True, plot total spectrum (sum over all components) from given ADF15 file.

• plot_all_lines (bool) – If True, plot all individual lines, rather than just the profiles due to different atomic processes. If more than 50 lines are included, a down-selection is automatically made to avoid excessive memory consumption.

• no_leg (bool) – If True, no plot legend is shown. Default is False, i.e. show legend.

• ax (matplotlib Axes instance) – Axes to plot on if plot=True. If left to None, a new figure is created.

Returns:

• wave_final_A (1D array) – Array of wavelengths in units of \(\r{A}\) on which the total spectrum is returned.

• spec_ion (1D array) – Spectrum from ionizing components of the input ADF15 file as a function of wave_final_A.

• spec_exc (1D array) – Spectrum from excitation components of the input ADF15 file as a function of wave_final_A.

• spec_rr (1D array) – Spectrum from radiative recombination components of the input ADF15 file as a function of wave_final_A.

• spec_dr (1D array) – Spectrum from dielectronic recombination components of the input ADF15 file as a function of wave_final_A.

• spec_cx (1D array) – Spectrum from charge exchange recombination components of the input ADF15 file as a function of wave_final_A.

• ax (matplotlib Axes instance) – Axes on which the plot is returned.

Notes

Including ionizing, excited and recombining charge states allows for a complete description of spectral lines that may derive from various atomic processes in a plasma.

Doppler broadening depends on the local ion temperature and mass of the emitting species. It is modeled here using

\[\theta(\nu) = \frac{1}{\sqrt{\pi}\Delta \nu_D} e^{-\left(\frac{\nu - \nu_0}{\Delta \nu_D}\right)^2}\]

with the Doppler profile half-width being

\[\Delta \nu_D = \frac{1}{\nu_0} \sqrt{\frac{2 T_i}{m}}\]

The Doppler shift dlam_A can be calculated from

\[\Delta \lambda_v = \lambda \cdot \left( 1 - \frac{v\cdot \cos(\alpha)}{c}\right)\]

where \(v\) is the plasma velocity and \(\alpha\) is the angle between the line-of-sight and the direction of plasma rotation.

Refs: S. Loch’s and C. Johnson’s PhD theses.

Minimal Working Example

Plot spectrum in one of the neutral H ADF15 files >>> filepath = aurora.get_adas_file_loc(‘pec96#h_pju#h0.dat’, filetype=’adf15’) >>> trs = aurora.read_adf15(filepath) >>> Te_eV = 80.; ne_cm3 = 1e14 # local ne [\(cm^{-3}\)]and Te [\(eV\)] >>> out = aurora.get_local_spectrum(filepath, ‘H’, ne_cm3, Te_eV) # defaults to excitation-only

aurora.radiation.get_main_ion_dens(ne_cm3, ions, rhop_plot=None)

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:

ni_cm3 – Estimated main ion density in \(cm^{-3}\).

Return type:

array (time,space)

aurora.radiation.get_photon_emissivity(adf15, lam, ne_cm3, Te_eV, imp_density, n0_cm3=0, meta_ind=None)

Evaluate PEC coefficients and return all components of intensity saved in an ADF15 file in units of \(ph/s\).

Parameters:

• adf15 (PandasFrame) – ADF15 PEC coefficients

• ne_cm3 (array (nr)) – Profile of electron density, in units of \(cm^{-3}\).

• Te_eV (array (nr)) – Profile of electron temperature, in units of \(eV\).

• imp_density (list or 2D array (nstate, nr)) – Density of all impurity states. These which are not needed can be set to None

• n0_cm3 (array, optional) – Local density of atomic neutral hydrogen isotopes. This is only used if the provided ADF15 file contains charge exchange contributions.

• meta_ind (list (nstate), optional) – indexes of metastates corresponding to imp_density

Returns:

intensity – arrays of line emissivity in ph/s units

Return type:

dict of arrays (nr)

aurora.radiation.parse_adf15_configs(path)

Parse ADF15 file to identify electron configurations and energies used to produce Photon Emissivity Coefficients (PECs) in the same file.

Parameters:

path (str) – Path to adf15 file to read.

Returns:

DataFrame containing information about all electron configurations used for PEC calculations.

Return type:

pandas.DataFrame

Notes

The format of ADAS ADF15 files has varied over time; consequently, it is surprisingly difficult to figure out parsing methods that work for all files. If this function fails on one of your files, please report this via a Github issue and one of the Aurora core developers will help adapting the parser.

aurora.radiation.parse_adf15_spec(lines, num_lines)

Parse full description of provided rates from an ADF15 file.

Parameters:

• lines (list or None) – List of lines read from ADF15 files.

• num_lines (int) – Number of spectral lines in the processed ADF15 file.

Returns:

Dictionary containing information about all loaded transitions.

Return type:

dict

aurora.radiation.plot_pec(transition, ax=None, plot_3d=False)

Plot a single Photon Emissivity Coefficient dependency on electron density and temperature.

Parameters:

• transition (pandas array with the PEC data) –

• ax (matplotlib axes instance) – If not None, plot on this set of axes.

• plot_3d (bool) – Display PEC data as 3D plots rather than 2D ones.

Examples

>>> filename = 'pec96#h_pju#h0.dat'
Fetch file automatically, locally, from AURORA_ADAS_DIR, or directly from the web:
>>> path = aurora.get_adas_file_loc(filename, filetype='adf15')
Load all transitions provided in the chosen ADF15 file:
>>> trs = aurora.read_adf15(path)
Select the Lyman-alpha transition and plot its rates:
>>> tr = trs[(trs['lambda [A]']==1215.2) & (trs['type']=='excit')]
>>> aurora.plot_pec(tr)
Alternatively, to select a line using the ADAS "ISEL" indices:
>>> aurora.plot_pec(trs.loc[(trs['isel']==1)])
or
>>> aurora.plot_pec(trs.iloc[0])

aurora.radiation.radiation_model(imp, rhop, ne_cm3, Te_eV, geqdsk, adas_files_sub={}, n0_cm3=None, Ti_eV=None, nz_cm3=None, frac=None, plot=False)

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

• geqdsk (dict, optional) – EFIT gfile as returned after postprocessing by the omfit_classes.omfit_eqdsk package (OMFITgeqdsk class).

• 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:

res – Dictionary containing results of radiation model.

Return type:

dict

aurora.radiation.read_adf15(path, order=1)

Read photon emissivity coefficients (PECs) from an ADAS ADF15 file.

Returns a pandas.DataFrame object describing all transitions available within the chosen ADF15 file.

PECs are provided in the form of an interpolant that will evaluate the log10 of the PEC at a desired electron density and temperature. The power-10 exponentiation of this PEC has units of \(photons \cdot cm^3/s\).

Units for interpolation: \(cm^{-3}\) for density; \(eV\) for temperature.

Parameters:

• path (str) – Path to adf15 file to read.

• order (int) – Parameter to control the order of interpolation. Default is 1 (linear interpolation).

Returns:

Parsed data from the given ADF15 file, including an interpolation function for the log-10 of the PEC of each spectral line. See the examples below for advice on using this form of output for plotting and further analysis.

Return type:

pandas.DataFrame

Examples

To plot the Lyman-alpha photon emissivity coefficients for H (or its isotopes), one can use:

>>> import aurora
>>> filename = 'pec96#h_pju#h0.dat'
>>> # fetch file automatically, locally, from AURORA_ADAS_DIR, or directly from the web:
>>> path = aurora.get_adas_file_loc(filename, filetype='adf15')
>>> # load all transitions provided in the chosen ADF15 file:
>>> trs = aurora.read_adf15(path)
>>> # select the excitation-driven component of the Lyman-alpha transition:
>>> tr = trs[(trs['lambda [A]']==1215.2) & (trs['type']=='excit')]
>>> # now plot the rates:
>>> aurora.plot_pec(tr)

Note that excitation-, recombination- and charge-exchange driven components can be fetched in the same way by specifying the type, e.g. using >>> trs[‘type’]==’excit’ # ‘excit’, ‘recom’ or ‘chexc’

Since the output of aurora.radiation.read_adf15 is a pandas DataFrame, its contents can be indexed in a variety of ways, e.g. one can select a line based on the ADAS “ISEL” indices and plot it using >>> aurora.plot_pec(trs.loc[(trs[‘isel’]==1)]) or simply >>> aurora.plot_pec(trs.iloc[0]) Spectral lines are ordered by ISEL numbers -1 (Python indexing!), so using the iloc method of a pandas DataFrame allows effective indexing by ISEL numbers.

The log-10 of the PEC interpolation function can be used as >>> tr[‘log10 PEC fun’].iloc[0].ev(np.log10(ne_cm3), np.log10(Te_eV)) where the interpolant was evaluated (via the ev method) at specific points of n_e (units of [\(cm^{-3}\)]) and T_e (units of \(eV\)). Note that the log-10 of n_e and T_e is needed, not n_e and T_e themselves!

Metastable-resolved files are automatically identified based on the file nomenclature and parsed accordingly, e.g.:

>>> filename = 'pec96#he_pjr#he0.dat'
>>> path = aurora.get_adas_file_loc(filename, filetype='adf15')
>>> trs = aurora.read_adf15(path)
>>> aurora.plot_pec(trs[trs['lambda [A]']==584.4])

Notes

This function expects the format of PEC files produced via the ADAS adas810 or adas218 routines.

aurora.radiation.sync_rad(B_T, ne_cm3, Te_eV, r_min, R_maj)

Calculate synchrotron radiation following Trubnikov’s formula [1]_. We make use of a simplified formulation as given by Zohm [2].

Parameters:

• B_T (float or 1D array) – Magnetic field amplitude [T].

• ne_cm3 (float or 1D array) – Electron density [\(cm^{-3}\)]

• Te_eV (float or 1D array) – Electron temperature [\(eV\)]

• r_min (float) – Minor radius [m].

• R_maj (float) – Major radius [m].

• Returns –

• array – Rate of synchrotron radiation [\(W/cm^3\)]

References

[1]

Trubnikov, JETP Lett. 16 25 (1972)

[2]

Zohm et al., Journal of Fusion Energy (2019) 38:3-10

aurora.grids_utils module

Methods to create radial and time grids for aurora simulations.

exception aurora.grids_utils.MissingAuroraBuild

Bases: UserWarning

aurora.grids_utils.create_aurora_time_grid(timing, plot=False)

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:

• time (array) – Computational time grid corresponding to timing input.

• save (array) – 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.

aurora.grids_utils.create_radial_grid(namelist, plot=False)

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:

• rvol_grid (array) – Volume-normalized grid used for Aurora simulations.

• pro (array) – Normalized first derivatives of the radial grid, defined as pro = (drho/dr)/(2 d_rho) = rho’/(2 d_rho)

• qpr (array) – Normalized second derivatives of the radial grid, defined as qpr = (d^2 rho/dr^2)/(2 d_rho) = rho’’/(2 d_rho)

• prox_param (float) – Grid parameter used for perpendicular loss rate at the last radial grid point.

aurora.grids_utils.create_radial_grid_fortran(namelist, plot=False)

This interfaces the package subroutine to create the radial grid exactly as STRAHL does it. Refer to the STRAHL manual for details.

aurora.grids_utils.create_time_grid(timing=None, plot=False)

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:

• time (array) – Computational time grid corresponding to :param:timing input.

• save (array) – 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.

aurora.grids_utils.create_time_grid_new(timing, verbose=False, plot=False)

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:

• t_vals (array) – Times in the time base [s]

• i_save (array) – Array of 0,1 values indicating at which times internal arrays should be stored/returned.

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

aurora.grids_utils.estimate_boundary_distance(shot, device, time_ms)

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:

• bound_sep (float) – Estimate for the distance between the wall boundary and the separatrix [cm]

• lim_sep (float) – 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.

aurora.grids_utils.estimate_clen(geqdsk)

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 omfit_classes.omfit_eqdsk.

Returns:

• clen_divertor (float) – Estimate of the connection length to the divertor [m]

• clen_limiter (float) – Estimate of the connection length to the limiter [m]

aurora.grids_utils.get_HFS_LFS(geqdsk, rho_pol=None)

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_rvol_mapping() for an application.

Parameters:

• geqdsk (dict) – Dictionary containing the g-EQDSK file as processed by the omfit_classes.omfit_eqdsk.

• 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:

• Rhfs (array) – Major radius [m] on the HFS

• Rlfs (array) – Major radius [m] on the LFS

aurora.grids_utils.get_rhopol_rvol_mapping(geqdsk, rho_pol=None)

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 omfit_classes.omfit_eqdsk.

• 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:

• rho_pol (array) – Sqrt of normalized poloidal flux grid

• rvol (array) – Mapping of rho_pol to a radial grid defined in terms of normalized flux surface volume.

aurora.grids_utils.vol_int(var, rvol_grid, pro_grid, Raxis_cm, rvol_max=None)

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:

• var (2D+ array (time, ..., radius)) – Data array for which a volume integral should be evaluated. The last dimension must be radial, other dimensions are arbitrary.

• rvol_grid (1D array) – Volume-normalized radial grid.

• pro_grid – Normalized first derivative of the radial grid, see create_radial_grid().

• Raxis_cm (float) – Major radius on axis [cm]

• rvol_max (float) – Maximum volume-normalized radius for integral. If not provided, integrate over the entire simulation radial grid.

Returns:

var_volint – Time history of volume integrated variable

Return type:

array (nt,)

aurora.coords module

aurora.coords.get_rhop_RZ(R, Z, geqdsk)

Find rhop at every R,Z [m] based on the equilibrium in the geqdsk dictionary.

aurora.coords.rV_vol_average(quant, r_V)

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 :math:`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.

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.

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)

Transform from one radial coordinate to another. Note that this coordinate conversion is only strictly valid inside of the LCFS. A number of common coordinate nomenclatures are accepted, but it is recommended to use one of the coordinate names indicated in the input descriptions below.

Parameters:

• x (array or float) – input x coordinate

• name_in (str) – input x coordinate name (‘rhon’,’rvol’,’rhop’,’rhov’,’Rmid’,’rmid’,’r/a’)

• name_out (str) – input x coordinate (‘rhon’,’psin’,’rvol’, ‘rhop’,’rhov’,’Rmid’,’rmid’,’r/a’)

• geqdsk (dict) – gEQDSK dictionary, as obtained from the omfit-eqdsk package.

Returns:

Conversion of x input for the requested radial grid coordinate.

Return type:

array

aurora.coords.rhoTheta2RZ(geqdsk, rho, theta, coord_in='rhop', n_line=201)

Convert values of rho,theta into R,Z coordinates from a geqdsk dictionary.

Parameters:

• geqdsk (output of the omfit_classes.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_classes.omfit_eqdsk. In this case, device, shot and time to fetch the equilibrium

• rho (np.ndarray) – Values of normalized radial coordinate to consider.

• theta (np.ndarray) – Values of poloidal angle coordinate to consider.

• coord_in (str) – Label describing the nature of the radial coordinate in use.

• n_line (int) – Number of points to discretize flux surface.

• Results –

• ------- –

• R (np.array, (ntheta, nrho)) – Values of the major radius along flux surfaces.

• Z (np.array, (ntheta, nrho)) – Values of the vertical coordinate along flux surfaces.

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

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_classes.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_classes.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_classes.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_classes.omfit_eqdsk dictionary is also returned

Returns:

• quant_vol_avg (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.

• geqdsk (dict) – Only returned if return_geqdsk=True.

aurora.source_utils module

Methods related to impurity source functions.

aurora.source_utils.get_radial_source(namelist, rvol_grid, pro_grid, S_rates, Ti_eV=None)

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.

• Ti_eV (array, optional (nt,nr)) – 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:

source_rad_prof – Radial profile of the impurity neutral source for each time step.

Return type:

array (nr,nt)

aurora.source_utils.get_source_time_history(namelist, Raxis_cm, time)

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

Parameters:

• namelist (dict) – Aurora namelist dictionary. The field namelist[‘source_type’] specifies how the source function is being specified – see the notes below.

• 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:

source_time_history – The source time history on the input time base.

Return type:

array (nt,)

Notes

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

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

#. namelist[‘source_type’] == ‘interp’: the time history for the source is provided by the user within the ‘explicit_source_time’ and ‘explicit_source_vals’ fields of the namelist dictionary and this data is simply interpolated.

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

#. 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’]

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

aurora.source_utils.lbo_source_function(t_start, t_rise, t_fall, n_particles=1.0, time_vec=None)

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:

• time_vec (array) – Times for the source function of each given impurity

• source (array) – Time history of the synthetized source function.

aurora.source_utils.read_source(filename)

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

Parameters:

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

Returns:

• t (array of float, (n,)) – The timebase (in seconds).

• s (array of float, (n,)) – The source function (#/s).

aurora.source_utils.write_source(t, s, shot, imp='Ca')

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:

contents – Content of the source file written to {imp}flx{shot}.dat

Return type:

str

aurora.plot_tools module

class aurora.plot_tools.DraggableColorbar(cbar, mapimage)

Bases: object

Create a draggable colorbar for matplotlib plots to enable quick changes in color scale.

Example::

fig,ax = plt.subplots()
cntr = ax.contourf(R, Z, vals)
cbar = plt.colorbar(cntr, format='%.3g', ax=ax)
cbar = DraggableColorbar(cbar,cntr)
cbar.connect()

connect()

Matplotlib connection for button and key pressing, release, and motion.

disconnect()

key_press(event)

Key pressing event

on_motion(event)

Move if the mouse is over the colorbar.

on_press(event)

Button pressing; check if mouse is over colorbar.

on_release(event)

Upon release, reset press data

aurora.plot_tools.get_color_cycle(num=None, map='plasma')

Get an iterable to select different colors in a loop. Efficiently splits a chosen color map into as many (num) parts as needed.

aurora.plot_tools.get_line_cycle()

aurora.plot_tools.get_ls_cycle()

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

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

aurora.default_nml.load_default_namelist()

Load default namelist. Users should modify and complement this for a successful forward-model 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)

aurora.interp.exppol1(params, d, rLCFS, r)

aurora.interp.funct(params, rLCFS, r)

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)

Function ‘funct2’ in STRAHL manual.

aurora.interp.interp(x, y, rLCFS, r)

Function ‘interp’ used in STRAHL for D and V profiles.

aurora.interp.interp_quad(x, y, d, rLCFS, r)

Function ‘interp’ used for kinetic profiles.

aurora.interp.interpa_quad(x, y, rLCFS, r)

Function ‘interpa’ used for kinetic profiles

aurora.interp.ratfun(params, d, rLCFS, r)

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)

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.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.Lya_to_neut_dens(emiss_prof, ne, Te, ni=None, plot=True, rhop=None, rates_source='adas', axs=None)

Estimate ground state neutral density from measured emissivity profiles. This ignores possible molecular dynamics and effects that may be captured via forward modeling of neutral transport.

Parameters:

• emiss_prof (1D array) – Emissivity profile, units of \(W/cm\)

• ne (1D array) – Electron density, units of \(cm^{-3}\)

• Te (1D array) – Electron temperature, units of \(eV\)

• ni (1D array) – Main ion (H/D/T) density, units of \(cm^{-3}\). If left to None, this is internally set to ni=ne.

• plot (bool) – If True, plot some of the key density profiles.

• rhop (1D array) – Sqrt of normalized poloidal flux radial coordinate. Used only for plotting.

• rates_source (str) – Source of atomic rates. Possible choices are ‘adas’ or ‘colrad’

• axs (Axes instance) – If given, plot on these axes.

Returns:

N1 – Radial profile of estimated ground state atomic neutral density on the same grid as the input arrays. Units of \(cm^{-3}\).

Return type:

1D array

Examples

>>> N2_colrad,axs = Lya_to_neut_dens_basic(emiss_prof, ne, Te, ni,
>>>                                  plot=True, rhop=rhop, rates_source='colrad')

>>> N2_adas,axs = Lya_to_neut_dens_basic(emiss_prof, ne, Te, ni, plot=True, rhop=rhop,
>>>                       rates_source='adas',axs=axs)

aurora.neutrals.download_ehr5_file()

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)

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()

plot(field='Seff', fig=None, axes=None)

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

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(\frac{N_m^i}{N_1} \right) N_m + \left(\frac{N_m^{ii}}{n_i} \right) 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.

Parameters:

• m (int) – Principal quantum number of excited state of interest. 2<m<10

• N1 (float, list or 1D-array [\(cm^{-3}\)]) – Density of ions in the ground state. This must have the same shape as ni!

• ni (float, list or 1D-array [\(cm^{-3}\)]) – Density of ions corresponding to the atom under consideration. This must have the same shape as N1!

• ne (float, list or 1D-array [\(cm^{-3}\)]) – Electron density to evaluate atomic rates at.

• Te (float, list or 1D-array [\(eV\)]) – Electron temperature to evaluate atomic rates at.

• rad_prof (list, 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_label (str) – When rad_prof is not None, this is the label for the radial coordinate.

• plot (bool) – Display the excited state ratio

Returns:

Nm – Density of electrons in excited state n [\(cm^{-3}\)]

Return type:

array of shape [len(ni)=len(N1),len(ne),len(Te)]

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)

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.

• 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:

Ns – 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:

list of arrays

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)

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_keV, E_beam, m_bckg, m_beam, n_level)

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_keV (float, 1D or 2D array) – Target temperature [keV]. Results will be computed for each Ti_keV 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:

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

Return type:

float, 1D or 2D array

aurora.nbi_neutrals.get_NBI_imp_cxr_q(neut_fsa, q, rhop_kp, times_kp, Ti_eV, ne_cm3, include_fast=True, include_halo=True, debug_plots=False)

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 and ne may be time-dependent, with a time base given by times_kp, the FSA neutrals are expected to be time-independent. Hence, the resulting CXR rates will only have time dependence that reflects changes in Ti and ne, 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_kp (array-like) – Sqrt of poloidal flux radial coordinate for Ti profiles.

• times_kp (array-like) – Time base on which Ti_eV is given [s].

• Ti_eV (array-like) – Ion temperature profile on the rhop_kp, times_kp bases, in units of \(eV\).

• ne_cm3 (array-like) – Electron density profile on the rhop_kp, times_kp bases, in units of \(cm^{-3}\).

• 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:

rates – 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:

dict

Notes

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

aurora.nbi_neutrals.get_neutrals_fsa(neutrals, geqdsk, debug_plots=True)

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 (dictionary output of omfit_classes.omfit_eqdsk.OMFITgeqdsk class) – gEQDSK post-processed dictionary, as given by omfit_classes.omfit_eqdsk.

• debug_plots (bool, optional) – If True, various plots are displayed.

Returns:

neut_fsa – 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:

dict

aurora.nbi_neutrals.rotation_matrix(alpha, beta, gamma)

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_keV, m_i, m_n, n_level)

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_keV (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.

Returns:

rate – output reaction rate in [\(cm^3/s\)] units

Return type:

float or 1D array

Notes

This does not currently account for the effect of rotation! Doing so will require making the integration in this function 2-dimensional.

aurora.nbi_neutrals.uvw_xyz(u, v, w, origin, R)

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)

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.

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

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:

sigma – Cross section of selected process, in [\(cm^2\)] units.

Return type:

float

Notes

See comments in Janev & Smith 1993 for uncertainty estimates.

aurora.janev_smith_rates.js_sigma_cx_n1_q1(E)

Electron capture cross section for

\[H^{+} + H(1s) --> H + H^+\]

Notes

Section 2.3.1 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n1_q2(E)

Electron capture cross section for

\[He^{2+} + H(1s) --> He^+ + H^+\]

Notes

Section 3.3.1 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n1_q4(E)

Electron capture cross section for

\[Be^{4+} + H(1s) --> Be^{3+} + H^+\]

Notes

Section 4.3.1 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n1_q5(E)

Electron capture cross section for

\[B^{5+} + H(1s) --> B^{4+} + H^+\]

Notes

Section 4.3.2 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n1_q6(E)

Electron capture cross section for

\[C^{6+} + H(1s) --> C^{5+} + H^+\]

Notes

Section 4.3.3 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n1_q8(E)

Electron capture cross section for

\[O^{8+} + H(1s) --> O^{7+} + H^+\]

Notes

Section 4.3.4 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n1_qg8(E, q)

Electron capture cross section for

\[A^{q+} + H(1s) --> A^{(q-1)+} + H^+, q>8\]

Notes

Section 4.3.5, p.172, of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_n2_q2(E)

Electron capture cross section for

\[He^{2+} + H(n=2) --> He^+ + H^+\]

Notes

Section 3.3.2 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_ng1_q1(E, n1)

Electron capture cross section for

\[H^{+} + H(n) --> H + H^+ , n>1\]

Notes

Section 2.3.2 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_ng1_qg3(E, n1, q)

Electron capture cross section for

\[A^{q+} + H^*(n) --> A^{(q-1)+}+H^+ , q>3\]

Notes

Section 4.3.6, p.174, of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_cx_ng2_q2(E, n1)

Electron capture cross section for

\[He^{2+} + H*(n) --> He^+ + H^+ , n>2\]

Notes

Section 3.2.3 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.js_sigma_ioniz_n1_q8(E)

Ionization cross section for

\[O^{8+} + H(1s) --> O^{8+} + H^+ +e^-\]

Notes

Section 4.2.4 of Janev & Smith, NF 1993.

aurora.janev_smith_rates.plot_js_sigma(q=18)

Plot/check sensitivity of JS cross sections to beam energy. NB: cross section is taken to only depend on partially-screened ion charge

aurora.synth_diags module

aurora.synth_diags.centrifugal_asymmetry(rhop, Rlfs, omega, Zeff, A_imp, Z_imp, Te, Ti, main_ion_A=2, plot=False, nz=None, geqdsk=None)

Estimate impurity poloidal asymmetry effects from centrifugal forces.

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

\[n(r,\theta) = n_0(r) \times \exp\left(\lambda(\rho) (R(r,\theta)^2- R_0^2)\right)\]

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

Parameters:

• rhop (array (nr,)) – Sqrt of normalized poloidal flux grid.

• Rlfs (array (nr,)) – Major radius on the Low Field Side (LFS), at points corresponding to rhop values

• 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).

• A_imp (float) – Impurity ion atomic mass number (e.g. 40 for Ca)

• Z_imp (array (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.

• Te (array (nr,nt)) – Electron temperature (eV)

• Ti (array (nr, nt)) – Background ion temperature (eV)

• main_ion_A (int, optional) – Background ion atomic mass number. Default is 2 for D.

• plot (bool) – If True, plot asymmetry factor \(\lambda\) vs. radius and show the predicted 2D impurity density distribution at the last time point.

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

• geqdsk (dict) – Dictionary containing the omfit_classes.omfit_eqdsk reading of the EFIT g-file.

Returns:

CF_lam – Asymmetry factor, defined as \(\lambda\) in the expression above.

Return type:

array (nr,)

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

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.

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

aurora.kn1d module

Aurora functionality to set up and run KN1D to extract atomic and neutral background densities at the edge.

KN1D is a 1D kinetic neutral code originally developed by B.LaBombard (MIT). For information, refer to the KN1D manual .

Note that this Aurora module is merely a wrapper of KN1D. Users require an IDL license on the computer where this module is called in order to be able to run KN1D. The IDL (and Fortran) code themselves are automatically downloaded and compiled by this module.

aurora.kn1d.plot_emiss(res, check_collrad=True)

Plot profiles of Ly-a and D-alpha emissivity from the KN1D output. KN1D internally computes Ly-a and D-alpha emission using the Johnson-Hinnov coefficients; here we check the result of that calculation and compare it to the prediction from atomic data from the COLLRAD collisional-radiative model included in DEGAS2.

Parameters:

• res (dict) – Output dictionary from function run_kn1d().

• check_collrad (bool) – If True, compare KN1D prediction of Ly-a and D-a emission using Johnson-Hinnov rates using rates from COLLRAD.

aurora.kn1d.plot_exc_states(res)

Plot excited state fractions of atomic neutral density from a KN1D run.

Parameters:

res (dict) – Output dictionary from function run_kn1d().

aurora.kn1d.plot_input_kin_prof(rmid_to_wall_m, ne_cm3, Te_eV, Ti_eV, innermost_rmid_cm, bound_sep_cm, lim_sep_cm)

Plot extent of kinetic profiles entering KN1D calculation

aurora.kn1d.plot_overview(res)

Plot an overview of a KN1D run, showing both kinetic profile inputs and a small selection of the outputs.

Parameters:

res (dict) – Output dictionary from function run_kn1d().

aurora.kn1d.plot_transport(res)

Make a simple set of plots of gradient scale lengths and effective diffusion coefficients from the KN1D output.

Parameters:

res (dict) – Output dictionary from function run_kn1d().

aurora.kn1d.run_kn1d(rhop, ne_cm3, Te_eV, Ti_eV, geqdsk, p_H2_mTorr, clen_divertor_cm, clen_limiter_cm, bound_sep_cm, lim_sep_cm, innermost_rmid_cm=5.0, mu=2.0, pipe_diag_cm=0.0, vx=0.0, collisions={}, kin_prof_exp_decay_SOL=False, kin_prof_exp_decay_LS=False, ne_decay_len_cm=[1.0, 1.0], Te_decay_len_cm=[1.0, 1.0], Ti_decay_len_cm=[1.0, 1.0], ne_min_cm3=1000000000000.0, Te_min_eV=1.0, Ti_min_eV=1.0, plot_kin_profs=False)

Run KN1D for the given parameters. Refer to the KN1D manual for details.

Depending on the provided options, kinetic profiles are extended beyond the Last Closed Flux Surface (LCxFS) and the Limiter Shadow (LS) via exponential decays with specified decay lengths. It is assumed that the given kinetic profiles extend from the core until at least the LCFS. All inputs are taken to be time-independent.

This function automatically checks if a KN1D repository is available; if it is not, it obtains it from the web and compiles the necessary code.

Note that an IDL license must be available. Aurora does not currently include a Python translation of KN1D – it only acts as a wrapper.

Parameters:

• rhop (1D array) – Sqrt of poloidal flux grid on which ne_cm3, Te_eV and Ti_eV are given.

• ne_cm3 (1D array) – Electron density on rhop grid [\(cm^{-3}\)].

• Te_eV (1D array) – Electron temperature on rhop grid [\(eV\)].

• Ti_eV (1D array) – Main ion temperature on rhop grid [\(eV\)].

• geqdsk (omfit_classes.omfit_eqdsk.OMFITgeqdsk class instance) – gEQDSK file as processed by the omfit_classes.omfit_eqdsk.OMFITgeqdsk class.

• p_H2_mTorr (float) – Pressure of molecular hydrogen-isotopes measured at the wall. This may be estimated from experimental pressure gauges. This variable effectively sets the amplitude of the neutral source at the edge. Units of \(mTorr\).

• clen_divertor_cm (float) – Connection length from the midplane to the divertor [\(cm\)].

• clen_limiter_cm (float) – Connection length from the midplane to the limiter [\(cm\)].

• bound_sep_cm (float) – Distance between the wall/boundary and the separatrix [\(cm\)].

• lim_sep_cm (float) – Distance between the limiter and the separatrix [\(cm\)].

• innermost_rmid_cm (float) – Distance from the wall to solve for. Default is 5 cm.

• mu (float) – Atomic mass number of simulated species. Default is 2.0 (D).

• pipe_diag_cm (float) – Diameter of the pipe through which H2 pressure is measured (see p_H2_mTorr variable). If left to 0, this diameter is effectively set to infinity. Default is 0.

• vx (float) – Radial velocity imposed on neutrals. This only has a weak effect usually. Default is 0 [\(cm/s\)].

• collisions (dict) – Collision terms flags. Set each to True or False. If any of the flags are not given, all collision terms are internally set to be active. Possible flags are ‘H2_H2_EL’,’H2_P_EL’,’H2_H_EL’,’H2_HP_CX’,’H_H_EL’,’H_P_CX’,’H_P_EL’,’Simple_CX’

• kin_prof_exp_decay_SOL (bool) – If True, kinetic profiles are set to exponentially decay over the SOL region.

• kin_prof_exp_decay_LS (bool) – If True, kinetic profiles are set to exponentially decay over the LS region.

• ne_decay_len_cm (list of 2 float) – Exponential decay lengths of electron density in the SOL and LS regions. Default is [1,1] \(cm\).

• Te_decay_len_cm (float) – Exponential decay lengths of electron temperature in the SOL and LS regions. Default is [1,1] \(cm\).

• Ti_decay_len_cm (float) – Exponential decay lengths of main ion temperature in the SOL and LS regions. Default is [1,1] \(cm\).

• ne_min_cm3 (float) – Minimum electron density across profile. Default is \(10^{12} cm^{-3}\).

• Te_min_eV (float) – Minimum electron temperaure across profile. Default is \(eV\).

• Ti_min_eV (float) – Minimum main ion temperaure across profile. Default is \(eV\).

• plot_kin_profs (bool) – If True, kinetic profiles input to KN1D are plotted.

Returns:

KN1D results and inputs, all collected into a dictionary. See example script for an illustration of using this.

Return type:

dict

Notes

For an example application, see the examples/aurora_kn1d.py script.

aurora.solps module

Aurora tools to read SOLPS results and extract atomic neutral density from EIRENE output. These enable examination of charge exchange recombination and radiation due to the interaction of heavy ions and thermal neutrals.

aurora.solps.apply_mask(triang, geqdsk, max_mask_len=0.4, mask_up=False, mask_down=False)

Function to apply basic masking to a matplolib triangulation. This type of masking is useful to avoid having triangulation edges going outside of the true simulation grid.

Parameters:

• triang (instance of matplotlib.tri.triangulation.Triangulation) – Matplotlib triangulation object for the (R,Z) grid.

• geqdsk (dict) – Dictionary containing gEQDSK file values as processed by omfit_classes.omfit_eqdsk.

• max_mask_len (float) – Maximum length [m] of segments within the triangulation. Segments longer than this value will not be plotted. This helps avoiding triangulation over regions where no data should be plotted, beyond the actual simulation grid.

• mask_up (bool) – If True, values in the upper vertical half of the mesh are masked. Default is False.

• mask_down (bool) – If True, values in the lower vertical half of the mesh are masked. Default is False.

Returns:

triang – Masked instance of the input matplotlib triangulation object.

Return type:

instance of matplotlib.tri.triangulation.Triangulation

aurora.solps.get_fort44_info(NDX, NDY, NATM, NMOL, NION, NSTRA, NCL, NPLS, NSTS, NLIM)

Collection of labels and dimensions for all fort.44 variables, as collected in the SOLPS-ITER 2020 manual.

aurora.solps.get_fort46_info(NTRII, NATM, NMOL, NION)

Collection of labels and dimensions for all fort.46 variables, as collected in the SOLPS-ITER 2020 manual.

aurora.solps.get_mdsmap()

Load dictionary allowing a mapping of chosen variables with SOLPS variable names on MDS+.

class aurora.solps.solps_case(*args, **kwargs)

Bases: object

Read SOLPS output, either from

1. output files on disk

2. an MDS+ tree

If arguments are provided with no keys, it is assumed that paths to b2fstate andd b2fgmtry are being provided (option #1). If keyword arguments are provided instead, we first check whether b2fstate_path and b2fgmtry_path are being given (option #1) or solps_id is given to load results from MDS+ (option #2). Other keyword arguments can be provided to specify the MDS+ server and tree to use (defaults are for AUG), and to provide a gEQDSK file.

Parameters:

• b2fstate_path (str (option #1)) – Path to SOLPS b2fstate output file.

• b2fgmtry_path (str (option #1)) – Path to SOLPS b2fgmtry output file.

• solps_id (str or int (option #2)) – Integer identifying the SOLPS run/shot to load from MDS+.

• server (str, optional (option #2)) – MDS+ server to load SOLPS results. Default is ‘solps-mdsplus.aug.ipp.mpg.de:8001’.

• tree (str (option #2)) – Name of MDS+ tree to load data from. Default is ‘solps’.

• geqdsk (str, optional (option #1 and #2)) – Path to the geqdsk to load from disk. Users may also directly provide an instance of the omfit_classes.omfit_geqdsk.OMFITgeqdsk class that contains the processed gEQDSK file. If not provided, the code tries to reconstruct a geqdsk based on the experiment name and time of analysis stored in the SOLPS output. If set to None, no geqdsk is loaded and functionality related to the geqdsk is not used.

Notes

The b2fstate and b2fgmtry parser expects the filenames to be simply named “b2fstate” and “b2fgmtry”, or else it may not recognize fields appropriately.

Minimal Working Example

import aurora so = aurora.solps_case(solps_id=141349) # load case 141349 from AUG MDS+; must have access to the server! so.plot2d_b2(so.data(‘ne’))

data(varname)

Fetch data either from files or MDS+ tree.

Parameters:

varname (str) – Name of SOLPS variable, e.g. ne,te,ti,dab2,etc.

eval_LOS(pnt1, pnt2, vals, npt=501, method='linear', plot=False, ax=None, label=None)

Evaluate the SOLPS output field along the line-of-sight (LOS) given by the segment going from point pnt1 to point pnt2 in 3D geometry.

Parameters:

• pnt1 (array (3,)) – Cartesian coordinates x,y,z of first extremum of LOS.

• pnt2 (array (3,)) – Cartesian coordinates x,y,z of second extremum of LOS.

• vals (array (self.data('ny'), self.data('nx'))) – Data array for a variable of interest.

• npt (int, optional) – Number of points to use for the path discretization.

• method ({'linear','nearest','cubic'}, optional) – Method of interpolation.

• plot (bool, optional) – If True, display variation of the field quantity as a function of the LOS path length from point pnt1.

• ax (matplotlib axes, optional) – Instance of figure axes to use for plotting. If None, a new figure is created.

• label (str, optional) – Text identifying the LOS under examination.

Returns:

array – Values of requested SOLPS output along the LOS

Return type:

(npt,)

find_gfile()

Identify the name of the gEQDSK file from the directory where b2fgmtry is also located.

find_xpoint()

Find location of the x-point in (R,Z) coordinates by using the fact that it is the only vertex shared by 8 cells.

get_3d_path(pnt1, pnt2, npt=501, plot=False, ax=None)

Given 2 points in 3D Cartesian coordinates, returns discretized R,Z coordinates along the segment that connects them.

Parameters:

• pnt1 (array (3,)) – Cartesian coordinates x,y,z of first path extremum (expected in [m]).

• pnt2 (array (3,)) – Cartesian coordinates x,y,z of second path extremum (expected in [m]).

• npt (int, optional) – Number of points to use for the path discretization.

• plot (bool, optional) – If True, display displaying result. Default is False.

• ax (matplotlib axes instance, optional) – If provided, draw 3d path on these axes. Default is to create a new figure and also draw the B2.5 grid polygons.

Returns:

• pathR (array (npt,)) – R [m] points along the path.

• pathZ (array (npt,)) – Z [m] points along the path.

• pathL (array (npt,)) – Length [m] discretization along the path.

get_b2_patches()

Get polygons describing B2 grid as a mp.collections.PatchCollection object.

get_poloidal_prof(vals, plot=False, label='', rhop=1.0, topology='LSN', ax=None)

Extract poloidal profile of a quantity “quant” from the SOLPS run. This function returns a profile of the specified quantity at the designated radial coordinate (rhop=1 by default) as a function of the poloidal angle, theta.

Note that double nulls (‘DN’) are not yet handled.

Parameters:

• vals (array (self.data('ny'), self.data('nx'))) – Data array for a variable of interest.

• plot (bool) – If True, plot poloidal profile.

• label (string) – Label for plot

• rhop (float) – Radial coordinate, in rho_p, at which to take poloidal surface. Default is 1 (LCFS)

• ax (matplotlib axes instance) – Axes on which poloidal profile should be plotted. If not given, a new set of axes is created. This is useful to possibly overplot profiles at different radii.

Returns:

• theta_rhop (1D array) – Poloidal grid measured in degrees from LFS midplane on which prof_rhop is given

• prof_rhop (1D array) – Mean poloidal profile at rhop on theta_rhop grid.

• prof_rhop_std (1D array) – Standard deviation of poloidal profile at rhop on the theta_rhop grid, based on variations within +/-dr_mm/2 millimeters from the surface at rhop.

get_radial_prof(vals, dz_mm=5, theta=0, label='', plot=False)

Extract radial profiles of a quantity “quant” from the SOLPS run. This function returns profiles on the low- (LFS) and high-field-side (HFS) midplane, as well as flux surface averaged (FSA) ones.

Parameters:

• vals (array (self.data('ny'), self.data('nx'))) – Data array for a variable of interest.

• dz_mm (float) – Vertical range [mm] over which quantity should be averaged near the midplane. Mean and standard deviation of profiles on the LFS and HFS will be returned based on variations of atomic neutral density within this vertical span. Note that this does not apply to the FSA calculation. Default is 5 mm.

• theta (float (0-360)) – Poloidal angle [degrees] at which to take radial profile, measured from 0 degrees at Outer Midplane. Default is 0 degrees

• label (string) – Optional string label for plot and legend. Default is empty (‘’)

• plot (bool) – If True, plot radial profiles.

Returns:

• rhop_fsa (1D array) – Sqrt of poloidal flux grid on which FSA profiles are given.

• prof_fsa (1D array) – FSA profile on rhop_fsa grid.

• rhop_LFS (1D array) – Sqrt of poloidal flux grid on which LFS profile (prof_LFS) is given.

• prof_LFS (1D array) – Mean LFS midpane profile on rhop_LFS grid.

• prof_LFS_std (1D array) – Standard deviation of LFS midplane profile on the rhop_LFS grid, based on variations within +/-dz_mm/2 millimeters from the midplane.

• rhop_HFS (1D array) – Sqrt of poloidal flux grid on which the midplane HFS profile (prof_HFS) is given.

• prof_HFS (1D array) – Mean HFS midplane profile on rhop_HFS grid.

• prof_HFS_std (1D array) – Standard deviation of HFS midplane profile on rhop_HFS grid, based on variations within +/-dz_mm/2 millimeters from the midplane.

load_data(fields=None, P_idxs=None, R_idxs=None, Rmin=None, Rmax=None, Pmin=None, Pmax=None)

Load SOLPS output for each of the needed quantities

Parameters:

• fields (list or array) – List of fields to fetch from SOLPS output. If left to None, by default uses [‘ne’,’Te’,’nn’,’Tn’,’nm’,’Tm’,’Ti’]

• P_idxs (list or array) – Poloidal indices to load.

• R_idxs (list or array) – Radial indices to load.

• Rmin (int or None.) – Minimum major radius index to load, if R_idxs is not given

• Rmax (int or None) – Maximum major radius index to load, if R_idxs is not given

• Pmin (int or None) – Minimum poloidal index to load, if P_idxs is not given

• Pmax (int or None) – Maximum poloidal index to load, if P_idxs is not given

load_eirene_mesh()

Load EIRENE nodes from the fort.33 file and triangulation from the fort.34 file

load_fort44()

Load result from one of the fort.44 file with EIRENE output on B2 grid.

Returns:

out – Dictionary containing a subdictionary with keys for each loaded field.

Return type:

dict

load_fort46()

Load result from fort.46 file with EIRENE output on EIRENE grid.

Returns:

out – Dictionary for each loaded file containing a subdictionary with keys for each loaded field from each file.

Return type:

dict

load_mesh_extra()

Load the mesh.extra file.

plot2d_b2(vals, ax=None, scale='log', label='', lb=None, ub=None, **kwargs)

Method to plot 2D fields on B2 grids. Colorbars are set to be manually adjustable, allowing variable image saturation.

Parameters:

• vals (array (self.data('ny'), self.data('nx'))) – Data array for a variable of interest.

• ax (matplotlib Axes instance) – Axes on which to plot. If left to None, a new figure is created.

• scale (str) – Choice of ‘linear’,’log’ and ‘symlog’ for matplotlib.colors.

• label (str) – Label to set on the colorbar. No label by default.

• lb (float) – Lower bound for colorbar. If left to None, the minimum value in vals is used.

• ub (float) – Upper bound for colorbar. If left to None, the maximum value in vals is used.

• kwargs – Additional keyword arguments passed to the PatchCollection class.

plot2d_eirene(vals, ax=None, scale='log', label='', lb=None, ub=None, replace_zero=True, **kwargs)

Method to plot 2D fields from EIRENE.

Parameters:

• vals (array (self.triangles)) – Data array for an EIRENE variable of interest.

• ax (matplotlib Axes instance) – Axes on which to plot. If left to None, a new figure is created.

• scale (str) – Choice of ‘linear’,’log’ and ‘symlog’ for matplotlib.colors.

• label (str) – Label to set on the colorbar. No label by default.

• lb (float) – Lower bound for colorbar. If left to None, the minimum value in vals is used.

• ub (float) – Upper bound for colorbar. If left to None, the maximum value in vals is used.

• replace_zero (boolean) – If True (default), replace all zeros in ‘vals’ with minimum value in ‘vals’

• kwargs – Additional keyword arguments passed to the tripcolor function.

plot_radial_summary(ls='o-b')

Plot a summary of radial profiles (ne, Te, Ti, nn, nm), at the inner and outer targets, as well as the outer midplane.

Parameters:

ls (str) –

plot_solps_2d_overview()

Display 2D views of most important SOLPS outputs.

plot_wall_geometry()

Method to plot vessel wall segment geometry from wall_geometry field in fort.44 file

species_id()

Identify species included in SOLPS run, both for B2 and EIRENE quantities. This is only designed to work in the “full” data format.

Module contents