Tutorial

Assuming that you have Aurora already installed on your system, we’re now ready to move forward. Some basic Aurora functionality is demonstrated in the examples package directory, where users may find a number of useful scripts. Here, we go through some of the same examples and methods.

Running Aurora simulations

If Aurora is correctly installed, you should be able to do:

import aurora

and then load a default namelist for impurity transport forward modeling:

namelist = aurora.load_default_namelist()

Note that you can always look at where this function is defined in the package by using, e.g.:

aurora.load_default_namelist.__module__

Once you have loaded the default namelist, have a look at the namelist dictionary. It contains a number of parameters that are needed for Aurora runs. Some of them, like the name of the device, are only important if automatic fetching of the EFIT equilibrium through MDSplus is required, or else it can be ignored (leaving it to its default value). Most of the parameter names should be fairly self-descriptive, but a detailed description will be available soon. In the meantime, please refer to docstrings through the code documentation.

Next, read in a magnetic equilibrium. You can find an example from a C-Mod discharge in the examples directory:

geqdsk = omfit_eqdsk.OMFITgeqdsk('example.gfile')

The output geqdsk dictionary contains the contents of the EFIT geqdsk file, with additional processing done by the omfit_eqdsk package for flux surfaces. Only some of the dictionary fields are used; refer to the grids_utils methods for details. The geqdsk dictionary is used to create a mapping between the rhop grid (square root of normalized poloidal flux) and a rvol grid, defined by the normalized volume of each flux surface. Aurora, like STRAHL, runs its simulations on the rvol grid.

We next need to read in some kinetic profiles, for example from an input.gacode file (available in the examples directory):

inputgacode = omfit_gapy.OMFITgacode('example.input.gacode')

Other file formats (e.g. plasma statefiles, TRANSP outputs, etc.) may also be read with omfit_gapy or other OMFIT-distributed packages. It is however not important to Aurora how the users get kinetic profiles: all that matters is that they are stored in the namelist[‘kin_prof’] dictionary. To set up time-independent kinetic profiles we can use:

kp = namelist['kin_profs']
kp['Te']['rhop'] = kp['ne']['rhop'] = np.sqrt(inputgacode['polflux']/inputgacode['polflux'][-1])
kp['ne']['vals'] = inputgacode['ne']*1e13    # 1e19 m^-3 --> cm^-3
kp['Te']['vals'] = inputgacode['Te']*1e3     # keV --> eV

Note that both electron density (ne) and temperature (Te) must be saved on a rhop grid. This grid is internally used by Aurora to map to the rvol grid. Also note that, unless otherwise stated, Aurora inputs are always in CGS units, i.e. all spatial quantities are given in \(cm\)!! (the extra exclamation mark is there for a good reason…).

Next, we specify the ion species that we want to simulate. We can simply do:

imp = namelist['imp'] = 'Ar'

and Aurora will internally find ADAS data for that ion (assuming that this is one of the common ones for fusion modeling). The namelist also contains information on what kind of source of impurities we need to simulate; here we are going to select a constant source (starting at t=0) of \(10^{24}\) particles/second.:

namelist['source_type'] = 'const'
namelist['Phi0'] = 1e24

Time dependent time histories of the impurity source may however be given by selecting namelist[‘source_type’]=”step” (for a series of step functions), “synth_LBO” (for an analytic function resembling a laser-blow-off (LBO) time history) or “file” (to load a detailed function from a file). Refer to the get_source_time_history() method for more details.

Assuming that we’re happy with all the inputs in the namelist at this point (many more could be changed!), we can now go ahead and set up our Aurora simulation::

asim = aurora.aurora_sim(namelist, geqdsk=geqdsk)

The aurora_sim class creates a Python object with spatial and temporal grids, kinetic profiles, atomic rates and all other inputs to the forward model. Aurora uses a diffusive-convective model for particle fluxes, so we need to specify diffusion (D) and convection (V) coefficients next::

D_z = 1e4 * np.ones(len(asim.rvol_grid))  # cm^2/s
V_z = -2e2 * np.ones(len(asim.rvol_grid)) # cm/s

Here we have made use of the rvol_grid attribute of the asim object, whose name is self-explanatory. This grid has a 1-to-1 correspondence with asim.rhop_grid. In the lines above we have created flat profiles of \(D=10^4 cm^2/s\) and \(V=-2\times 10^2 cm/s\), defined on our simulation grids. D’s and V’s could in principle (and, very often, in practice) be defined with more dimensions to represent a time-dependence and also different values for different charge states. Unless specifed otherwise, Aurora assumes all points of the time grid (now stored in asim.time_grid) and all charge states to have the same D and V. See the run_aurora() method for details on how to speficy further dependencies.

At this point, we are ready to run an Aurora simulation, with:

out = asim.run_aurora(D_z, V_z)

Blazing fast! Depending on how many time and radial points you have requested (a few hundreds by default), how many charge states you are simulating, etc., a simulation could take as little as <50 ms, which is significantly faster than other code, as far as we know. If you add use_julia=True to the run_aurora() call the run will be even faster; wear your seatbelt!

You can easily check the quality of particle conservation in the various reservoirs by using:

reservoirs = asim.check_conservation()

which will show the results in full detail. The reservoirs output list contains information about how many particles are in the plasma, in the wall reservoir, in the pump, etc.. Refer to the run_aurora() docstring for details.

A plot is worth a thousand words, so let’s make one for the charge state densities (on a nice slider!):

aurora.slider_plot(asim.rvol_grid, asim.time_out, asim.rad['line_rad'].transpose(1,2,0),
                   xlabel=r'$r_V$ [cm]', ylabel='time [s]', zlabel='Total radiation [A.U.]',
                   labels=[str(i) for i in np.arange(0,nz.shape[1])],
                   plot_sum=True, x_line=asim.rvol_lcfs )

Use the slider to go over time, as you look at the distributions over radius of all the charge states. It would be really great if you could just save this type of time- and spatially-dependent visualization to a video-format, right? That couldn’t be easier, using the animate_aurora() function::

aurora.animate_aurora(asim.rhop_grid, asim.time_out, nz.transpose(1,0,2),
                      xlabel=r'$\rho_p$', ylabel='t={:.4f} [s]', zlabel=r'$n_z$ [A.U.]',
                      labels=[str(i) for i in np.arange(0,nz.shape[1])],
                      plot_sum=True, save_filename='aurora_anim')

After running this, a .mp4 file with the name “aurora_anim.mp4” will be saved locally.

Radiation predictions

Once a set of charge state densities has been obtained, it is simple to compute radiation terms in Aurora. For example, using the results from the Aurora run in Running Aurora simulations, one can then run:

asim.rad = aurora.compute_rad(imp, nz.transpose(2,1,0), asim.ne, asim.Te, prad_flag=True)

The documentation on compute_rad() gives details on input array dimensions and various flags that may be turned on. In the case above, we simply indicated the ion number (imp), and provided charge state densities (with dimensions of time, charge state and space), electron density and temperature (dimensions of time and space). We then explicitely indicated prad_flag=True, which means that unfiltered “effective” radiation terms (line radiation and continuum radiation) should be computed. Bremsstrahlung is also estimated using an interpolation formula that is independent of ADAS data and can be found in asim.rad[‘brems’]. However, note that bremsstrahlung is already included in asim.rad[‘cont_rad’], which also includes other terms including continuum recombination using ADAS data. It can be useful to compare the bremsstrahlung calculation in asim.rad[‘brems’] with asim.rad[‘cont_rad’], but we recommend that users rely on the full continuum prediction for total power estimations.

Other possible flags of the compute_rad() function include:

  1. sxr_flag: if True, compute line and continuum radiation in the SXR range using the ADAS “pls” and “prs” files. Bremsstrahlung is also separately computed using the ADAS “pbs” files.

  2. thermal_cx_rad_flag: if True, the code checks for inputs n0 (atomic H/D/T neutral density) and Ti (ion temperature) and computes line power due to charge transfer from thermal background neutrals and impurities.

  3. spectral_brem_flag: if True, use the ADAS “brs” files to compute bremsstrahlung at a wavelength specified by the chosen file.

All of the radiation flags are False by default.

ADAS files for all calculations are taken by default from the list of files indicated in adas_files_dict() function, but may be replaced by specifying the adas_files dictionary argument to compute_rad().

Results from compute_rad() are collected in a dictionary (named “rad” above and added as an attribute to the “asim” object, for convenience) with clear keys, described in the function documentation. To get a quick plot of the radiation profiles, e.g. for line radiation from all simulated charge states, one can do:

aurora.slider_plot(asim.rvol_grid, asim.time_out, asim.rad['line_rad'].transpose(1,2,0),
                   xlabel=r'$r_V$ [cm]', ylabel='time [s]', zlabel='Total radiation [A.U.]',
                   labels=[str(i) for i in np.arange(0,nz.shape[1])],
                   plot_sum=True, x_line=asim.rvol_lcfs)

Aurora’s radiation modeling capabilities may also be useful when assessing total power radiation for integrated modeling. The radiation_model() function allows one to easily obtain the most important radiation terms at a single time slice, both as power densities (units of \(MW/cm^{-3}\)) and absolute power (units of \(MW\)). To obtain the latter form, we need to integrate over flux surface volumes. We can use the geqdsk dictionary obtained via:

geqdsk = omfit_eqdsk.OMFITgeqdsk('example.gfile')

(or equivalent methods/files) to then extract flux surface volumes (units of \(m^3\)) at each value of rhop::

grhop = np.sqrt(geqdsk['fluxSurfaces']['geo']['psin'])
gvol = geqdsk['fluxSurfaces']['geo']['vol']

# interpolate on our grid
vol = interp1d(grhop, gvol)(rhop)

We can now pass the vol array to radiation_model(), together with the impurity atomic symbol (imp), the rhop grid array, electron density (ne_cm3) and temperature (Te_eV) and, optionally, also background neutral densities to include thermal charge exchange::

res = aurora.radiation_model(imp,rhop,ne_cm3,Te_eV, vol,
                             n0_cm3=None, frac=0.005, plot=True)

Here we specified the impurity densities as a simple fraction of the electron density profile, by specifying the frac argument. This is obviously a simplifying assumption, effectively stating that the total impurity density profile should have a maximum amplitude of frac (in the case above, set to 0.005) and a profile shape (corresponding t a profile of V/D) that is identical to the one of the \(n_e\) profile. This may be convenient for parameter scans in the design process of future devices, but is by no means a correct assumption. If we’d rather calculate the total radiated power from a specific set of impurity charge state profiles (e.g. from an Aurora simulation), we can do:

res = aurora.radiation_model(imp,rhop,ne_cm3,Te_eV, vol,
                             n0_cm3=None, nz_cm3=nz_cm3, plot=True)

where we specified the charge state densities (dimensions of space, charge state) at a single time. Since we specified plot=True, a number of useful radiation profiles should be displayed.

Of course, one can also estimate radiation from the main ions. To do this, we first want to estimate the main ion density, using:

ni_cm3 = aurora.get_main_ion_dens(ne_cm3, ions)

with ions being a dictionary of the form:

ions = {'C': nC_cm3, 'Ne': nNe_cm3}   # (time,charge state,space)

with a number of impurity charge state densities with dimensions of (time,charge state,space). The get_main_ion_dens() function subtracts each of these densities (times the Z of each charge state) from the electron density to obtain a main ion density estimate based on quasineutrality. Before we move forward, we need to add a neutral stage density for the main ion species, e.g. using:

niz_cm3 = np.vstack((n0_cm3[None,:],ni_cm3)).T

such that the niz_cm3 output is a 2D array of dimensions (charge state, radius).

To estimate main ion radiation we can now do:

res_mainion = aurora.radiation_model('H',rhop,ne_cm3,Te_eV, vol, nz_cm3 = niz_cm3, plot=True)

(Note that the atomic data does not discriminate between hydrogen isotopes) In the call above, the neutral density has been included in niz_cm3, but note that (1) there is no radiation due to charge exchange between deuterium neutrals and deuterium ions, since they are indistinguishable, and (2) we did not attempt to include the effect of charge exchange on deuterium fractional abundances because n0_cm3 (included in niz_cm3 already fully specifies fractional abundances for main ions).

Zeff contributions

Following an Aurora run, one may be interested in what is the contribution of the simulated impurity to the total effective charge of the plasma. The calc_Zeff() method allows one to quickly compute this by running:

asim.calc_Zeff()

This makes use of the electron density profiles (as a function of space and time), stored in the “asim” object, and keeps Zeff contributions separate for each charge state. They can of course be plotted with slider_plot()::

aurora.slider_plot(asim.rvol_grid, asim.time_out, asim.delta_Zeff.transpose(1,0,2),
                   xlabel=r'$r_V$ [cm]', ylabel='time [s]', zlabel=r'$\Delta$ $Z_{eff}$',
                   labels=[str(i) for i in np.arange(0,nz.shape[1])],
                   plot_sum=True,x_line=asim.rvol_lcfs)

Ionization equilibrium

It may be useful to compare and contrast the charge state distributions obtained from an Aurora run with the distributions predicted by pure ionization equilibium, i.e. by atomic physics only, with no trasport. To do this, we only need some kinetic profiles, which for this example we will load from the sample input.gacode file available in the “examples” directory::

import omfit_gapy
inputgacode = omfit_gapy.OMFITgacode('example.input.gacode')

Recall that Aurora generally uses CGS units, so we need to convert electron densities to \(cm^{-3}\) and electron temperatures to \(eV\):

rhop = np.sqrt(inputgacode['polflux']/inputgacode['polflux'][-1])
ne_vals = inputgacode['ne']*1e13 # 1e19 m^-3 --> cm^-3
Te_vals = inputgacode['Te']*1e3  # keV --> eV

Here we also defined a rhop grid from the poloidal flux values in the inputgacode dictionary. We can then use the get_atom_data() function to read atomic effective ionization (“scd”) and recombination (“acd”) from the default ADAS files listed in adas_files_dict(). In this example, we are going to focus on calcium ions::

atom_data = aurora.get_atom_data('Ca',['scd','acd'])

In ionization equilibrium, all ionization and recombination processes will be perfectly balanced. This condition corresponds to specific fractions of each charge state at some locations that we define using arrays of electron density and temperature. We can compute fractional abundances and plot results using:

logTe, fz, rates = aurora.get_frac_abundances(atom_data, ne_vals, Te_vals, rho=rhop, plot=True)

The get_frac_abundances() function returns the log-10 of the electron temperature on the same grid as the fractional abundances, given by the fz parameter (dimensions: space, charge state). This same function can be used to both compute radiation profiles of fractional abundances or to compute fractional abundances as a function of scanned parameters ne and/or Te. The inverse of the rates output correspond to the atomic relaxation time. An additional argument of ne_tau (units of \(m^{-3}\cdot s\)) can be used to approximately model the effect of transport on ionization balance.

Working with neutrals

Aurora includes a number of useful functions for neutral modeling, both from the edge of fusion devices (thermal neutrals) and from neutral beams (fast and halo neutrals).

For thermal neutrals, we make use of atomic data from the Collrad collisional-radiative model, part of the DEGAS2 code.

The erh5_file class allows one to parse the erh5.dat file of DEGAS-2 that contains useful information to assess excited state fractions of neutrals in specific kinetic backgrounds. If the erh5.dat file is not available already, Aurora will download it and store it locally within its distribution directory. The data in this file is used for example in the get_exc_state_ratio() function, which given a ground state density of neutrals (N1), some ion and electron densities (ni and ne) and electron temperature (Te), will compute the fraction of neutrals in the principal quantum number m. Keyword arguments can be passed to this function to plot the results. Note that kinetic inputs may be given as a scalar or as a 1D list/array. The plot_exc_ratios() function may also be useful to plot the excited state ratios.

Note that in order to find the photon emissivity coefficient of specific neutral lines, the read_adf15() function may be used. For example, to obtain interpolation functions for neutral H Lyman-alpha emissivity, one can use:

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

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

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

This will plot the Lyman-alpha photon emissivity coefficients (both the components due to excitation and recombination) as a function of temperature in eV. Some files (e.g. try pec96#c_pju#c2.dat) may also have charge exchange components.

Analysis routines to work with fast and halo neutrals are also provided in Aurora. Atomic rates for charge exchange of impurities with NBI neutrals are taken from Janev & Smith NF 1993 and can be obtained from js_sigma(), which wraps a number of functions for specific atomic processes. To compute charge exchange rates between NBI neutrals (fast or thermal) and any ions in the plasma, users need to provide a prediction of neutral densities, likely from an external code like FIDASIM.

Neutral densities for each fast ion population (full-,half- and third-energy), multiple halo generations and a few excited states are expected. Refer to the documentation of get_neutrals_fsa() to read about how to provide neutrals on a poloidal cross section so that they may be “flux-surface averaged”.

bt_rate_maxwell_average() shows how beam-thermal Maxwell-averaged rates can be obtained; tt_rate_maxwell_average() shows the equivalent for thermal-thermal Maxwell-averaged rates.

Finally, get_NBI_imp_cxr_q() shows how flux-surface-averaged charge exchnage recombination rates between an impurity ion of charge q with NBI neutrals (all populations, fast and thermal) can be computed for use in Aurora forward modeling. For more details, feel free to contact Francesco Sciortino (sciortino-at-psfc.mit.edu).