Processing SOLPS-ITER output

SOLPS-ITER has several built-in tools for output visualisation and processing, the most prominent being b2plot. Sometimes, though, it is necessary to access the SOLPS-ITER output using more flexible software, like Matlab or Python. I suppose the SOLPS GUI counts as a mean of post-processing as well. This section introduces the basics of all post-processing methods we know.

Workflow fragmentation and remote access

Ideally, you have everything in one place: your SOLPS-ITER installation, experimental data and post-processing routines. Unfortunately, this is often not the case. At some point of my career, I was running SOLPS-ITER at Marconi Gateway, processed the output at IPP Garching and compared it to experimental data at IPP Prague. Further complicating matters, SOLPS-ITER requires a special work environment, SOLPSpy requires Python 2, current COMPASS modules run Python 3 of various versions and Matlab is a chapter of its own.

If you run SOLPS-ITER in a different place than your data post-processing, the best way to access the code output is to mount the location of your SOLPS-ITER files to wherever you do the post-processing. If you for some reason can't mount (like on AUG workstations where mounting requires root privileges), the next option is to hard copy the entire run folder, for instance by the rsync command. See more in the Remote access document, particularly section Transferring SOLPS-ITER runs. For some post-processing methods, copying only a few select files will do. For solps-postproc, for instance, these are:

b2fstate
b2fplasmf
b2fgmtry  # copied from baserun
input.dat
b2mn.dat

B2plot

b2plot is the standard built-in tool for visualising SOLPS-ITER results, and boy, over decades of refinement it has become buff like Arnold Schwarzenegger. It can do some really amazing things, provided you can figure out how from the not-very-telling section in the manual. Here I give a couple of beginner tips. Further information can be found in Questions and answers: b2plot and in the manual appendix I: b2plot manual.

The basics

To use b2plot, first get inside your selected run folder and, in the SOLPS-ITER environment, type:

echo "te surf" | b2plot

This will plot the electron temperature (te) by filling the appropriate cell with colour (surf), like so:

It looks so weird because this is the "computation domain" - the grid rectangles resting one next to the other. To see the "physical domain", you have to add the keyword phys to your command chain:

echo "phys te surf" | b2plot

Voilà. Instead of the electron temperature, you may plot any other B2.5 quantity in the same manner: the electron density (ne), all of the ion densities (na), ion temperature (ti) and many more listed in the I.2 section of the manual.

Data stacks

As the manual appendix I: b2plot manual explains, the data which b2plot can display come in seven different stacks:

stack name	stack content
matrix stack	B2.5 results
integer stack	integer arguments for b2plot commands
real stack	real arguments for b2plot commands
string stack	string arguments for b2plot commands
EIRENE stack	EIRENE results
wall stack	SOLPS-ITER results which are defined along the wall (such as sputtering yields)
target stack	SOLPS-ITER results which are defined along the divertor targets (such as heat loads)

This is important because b2plot is very peculiar about its stacks and many of its error messages mention them. For example, if you run

`echo "phys te 40 fmax surf" | b2plot`

...b2plot will complain about Insufficient data on real stack and ignore the fmax function (setting the colour bar maximum). This is because fmax takes a real number argument (40.0), not an integer (40). This is stated in the manual, but it's common error which can be fixed on-the-fly if you know what the error message means.

Another situation where you need to know about stacks is when plotting EIRENE output or making wall/target plots. Adapting the minimum example above, you may attempt to run

echo "phys edena surf" | b2plot

...to plot the EIRENE atom energy density edena (proportional to the atom pressure). However, this will result in the error message Insufficient data on matrix stack and no picture being drawn. The reason is that the surf command only plots data from the matrix stack (B2.5 results, defined on the rectangular B2.5 grid). To plot data from the EIRENE stack, defined on the triangular EIRENE grid, you need the esurf command instead.

In a similar vein, logf turns on the logarithmic scale for matrix data, while logw turns it on for wall data. zsel selects the species index for matrix data, while wzsel selects it for wall data. Thus the need for the beginner to know about the b2plot data stacks.

Passing arguments to b2plot commands

Always keep two things mind:

• Argument value precedes the command. For instance, b2plot will not understand zmin -0.3 but only -0.3 zmin or -3e-1 zmin.
• Some commands take integer arguments (3) and others float/real arguments (3.0). For instance, b2plot will not understand -1 zmin but only -1.0 zmin. Which format is appropriate for which command is given in the manual, section I.4 Available Functions and Operators in b2plot. For instance:
  • <number,I> font means the font command takes one integer argument, such as -25 font.
  • <label, S> glbl means the glbl command takes a string argument, such as 'Ugly carbon case' glbl.
  • <factor,R> lbsz means the lbsz command takes real number argument, such as 2.0 lbsz.

Composing b2plot commands and writing scripts

The longer a b2plot command is, the worse it is contained in the command line. Eventually, you may want to write a script, such as the following:

# global plot necessities
phys a4l 2 1 page vesl sep

# figure 1 wall plot: chemical sputtering
wall logw sputchem 2 wzsel 5.0 wwid

# figure 1 surf plot: C0 poloidal flux
fnax 2 zsel surf

# figure 2 wall plot: physical sputtering
wall logw sput 2 wzsel 5.0 wwid

# figure 2 surf plot: C0 poloidal flux
fnax 2 zsel surf

To run this script, located in the run directory and called sputtering.cmd:

echo '"sputtering.cmd" cmd' | b2plot

Let's break the script down. The first line, phys a4l 2 1 page vesl sep, defines the overall figure composition and common elements of its subfigures.

• phys: the individual subfigures will be plotted in the "physical domain" - that is, the axes are \(R\) and \(Z\) [m], and not \(x\) and \(y\) [cell number]
• a4l: the figure is A4 in format, in the landscape orientation
• 2 1 page: there are two subfigures in the figure, next to each other (2 columns, 1 row)
• vesl: in all the subfigures, the vacuum vessel will be plotted
• sep: in all the subfigures, the separatrix will be plotted

The second line, wall logw sputchem 2 wzsel 5.0 wwid, contains a wall plot (the lavender line around the B2.5 computational domain in the left subfigure):

• wall: a wall plot shall be drawn on top of a regular (here matrix) plot
• logw: set the colour bar scale to logarithmic
• sputchem: plot the chemical sputtering of all ion species (refer to Reactions of the Hexaaqua Ions with Carbonate Ions, The Hydronium Ion and [Ballauf 2020] for ion chemical sputtering)
• 2 wzsel: plot only ion species with index 2 (here C⁰)
• 5.0 wwid: set the line thickness to 5.0

The third line, fnax 2 zsel surf, contains a matrix plot underlying the wall plot:

• fnax: plot the poloidal particle flux of all ion species
• 2 zsel: plot only ion species with index 2 (here C⁰)
• surf: conclude drawing the subfigure by taking all the matrix stack data you have and make a surface plot out of them on the B2.5 grid

The fourth and fifth line are analogous to the second and third, only they concern the physical sputtering.

Notice that each subfigure has two colour bars: one for the wall plot (sputtering) and one for the matrix plot (C⁰ poloidal flux). Also notice that the wall plot shows "less than 10^-10" across the entire "wall" zone. This is a problem of the simulation set up (or a bug, depending on how smart you'd like b2plot to be). The sputchem and sput commands return the sputtering fluxes calculated by B2.5, while this is a coupled run where the sputtering is handled by EIRENE, so the B2.5 sputtering fluxes are zero.

From trial and error, it seems like b2plot favours a certain structure of its commands. To create more complicated scripts, try following the structure described above:

• global figure settings (size, orientation, axes labels, font size, plot area...)
• plotted quantity/variable (electron temperature, EIRENE molecules density...)
• commands pertaining to this particular variable/plot (axes extent, colour bar extent, colour palette, linear/logarithmic scale...)
• plot command (typically surf)

Note that some b2plot commands, such as mesh, have a built-in plot command and you needn't follow them up with surf.

echo "phys mesh" | b2plot

A few b2plot scripts

2D map of total carbon line radiation

a4p lbof logf phys vesl sep
b2ra pvol abs 2 8 sumz 1.0e3 fmin 1.0e7 fmax surf

2D map of all carbon species densities

phys a4l 4 2 page logf aspc vesl 1.5 lbsz
0.3 rmin 0.8 rmax -0.4 zmin 0.4 zmax
na 2 zsel 9e3 fmin 1.1e4 fmax surf
na 3 zsel 1e3 fmin 1e17 fmax surf
na 4 zsel 1e3 fmin 1e17 fmax surf
na 5 zsel 1e3 fmin 1e17 fmax surf
na 6 zsel 1e3 fmin 1e17 fmax surf
na 7 zsel 1e3 fmin 1e17 fmax surf
na 8 zsel 1e3 fmin 1e17 fmax surf

2D map of total neutral pressure (with B2.5 stacks)

Notice that since B2.5 stacks are used (rm*, surf etc.), the pressure is only given on the B2.5 grid.

# global plot necessities
phys a4p vesl sep

# multiply D atomic density and D atomic temperature to get atomic pressure
dab2 0 zsel tab2 0 zsel m*

# multiply molecular density and molecular temperature to get molecular pressure
dmb2 tmb2 m*

# add the pressures to get total neutral pressure
m+

# multiply by the electron charge to convert the units to Pa
1.6e-19 rm*

# set a larger label size
1.5 lbsz

# set colorbar extent
0. 2. 21 clvn

# finish the plot
surf

2D map of total neutral pressure (with EIRENE stacks)

Notice that since EIRENE stacks are used (rem*, esurf etc.), the pressure is given on the entire EIRENE grid.

# global plot necessities
phys a4p vesl sep

# add EIRENE neutral atom and molecule energy density to get the total neutral energy density in eV/m^3
edena edenm em+

# multiply by the electron charge to get the total neutral energy density in J/m^3
1.602e-19 rem*

# multiply by 2/3 to get the total neutral pressure in Pa
0.66 rem*

# set colorbar extent
0. 2. 21 eclvn

# finish the plot
esurf

B2.5 and EIRENE grids

This one's tricky because b2plot is stupid about choosing mesh colours. As far as I can tell, the mcol argument changes the colour of both the EIRENE grid (plotted by trig or trigmesh) and the B2.5 grid (plotted with mesh), so you can't plot them into one picture with different colours. My workaround is to plot the two grids into separate figures and overlay them in an external editor. The picture will look much better if you use the b2plot.ps file and a vector graphics editor (such as Inkscape) than if you convert the b2plot.ps file to a bitmap (see above) and use a raster graphics editor (such as GIMP).

# global plot necessities: first picture
phys a4p vesl sep

# plot the B2.5 + EIRENE grid in lavender
2 mcol trigmesh

# global plot necessities: second picture
phys a4p

# plot the B2.5 grid in blue
5 mcol mesh

Useful b2plot commands and functions

By adding words to the command chain, you can tweak the plot in many ways:

• adjust the font size for all text (default size is 1.0)

  echo "phys 2.0 lbsz te surf" | b2plot
  

• crop the plotted area (zmin and zmax in the z-direction, rmin and rmax in the r-direction)

  echo "phys -0.3 zmin 0.3 zmax te surf" | b2plot
  

  Note: In an effort to produce a pretty picture, b2plot engages in heavy optimisation calculating the exact axis bounds when the aspect is set to 1 (meaning 1 m on the X axis has the same length as 1 m on the Y axis). It tries to make the axis labels pretty, fill out the figure as best it can and preserve the aspect ratio. As a result, it can be very stubborn about the axis extent it has chosen and convincing it otherwise may be hard. If careful optimisation fails, I'd recommend using the toggle aspc, which toggles off the unity aspect ratio, and setting the figure size and axes extents so that the aspect ratio is approximately unity. No one's going to put a ruler to the picture anyway.

• plot the separatrix

  echo "phys sep te surf" | b2plot
  

• plot the vacuum vessel outline

  echo "phys vesl te surf" | b2plot
  

• set the axis labels (introduced in the 44th SOLPS-ITER User Forum)

  echo "phys te 'R [m]' xlab 'Z [m]' ylab surf" | b2plot
  

• set the plot title

  echo "phys te 'Electron temperature' glbl surf" | b2plot
  

• change the colourmap scale to logarithmic (useful for seeing more detail, but take care that your data doesn't take on negative values)

  echo "phys na logf surf" | b2plot
  

• set the colour bar limits and number of levels

  echo "phys pdena 0.0 1e19 32 eclvn esurf" | b2plot
  

• set the colour bar pallete

  echo "phys pdena 'blue_to_red_64.pal' cltb esurf" | b2plot
  

• plot something against the parallel coordinate

  echo "sx xpar 18 f.x" | b2plot
  

  This, in particular, plots the first SOL flux tube poloidal cross-section. Beware that grid indexing in b2plot begins from -1, not from 0!

Plotting EIRENE output

EIRENE data uses a whole different set of plot commands, usually with "e" prepended. To plot, for instance, the total atom density:

echo "phys edena esurf" | b2plot

Note that the command surf was replaced by esurf, that is, b2plot now uses the EIRENE grid.

A more complicated command to plot the total neutral pressure (in D-only cases, with impurities you'd also have to select ion species numbers):

echo "phys sep vesl logf edena edenm em+ 0.66 rem* esurf" | b2plot

Beside the plot commands specified above (phys, sep, vesl, logf and esurf), notice the operation on EIRENE stacks:

• Sum two EIRENE arrays: edena edenm em+. Notice the order: first the two arrays, then the EIRENE plus operation.
• Multiply an EIRENE array with a real number: edena edemn em+ 0.66 rem*. Again, notice the order "argument argument operation", and also that the number includes the decimal point. If the arrays were multiplied by 10, 10.0 would work while 10 would prompt the error message Insufficient data on real stack, since the code expected a real number and got an integer.

wlld

wlld is a very useful command, which extracts profiles of various quantities (most notoriously, heat fluxes broken down into individual components) along the divertor targets and the outer midplane. To create the text files containing these profiles, write:

echo "1100 wlld" | b2plot

This creates the files as ld tg *.dat, fp tg *.dat in the folder b2pl.exe.dir/. See the manual for more information.

Saving the plot results

The b2plot output is automatically saved under the name b2plot.ps (ps = PostScript, a vector graphic format). It's located either directly in the run/ folder or in its subfolder b2pl.exe.dir. I am not sure what the difference is, but surely it is in one of them. Every time you call b2plot again, the files are overwritten, so if you want to keep something, you have to rename it or make a copy of it.

While the .ps format is compatible with many devices, you may find raster formats (.png etc.) more convenient. In Linux distributions, the conversion may be done using the ImageMagick programme:

convert -density 150 b2plot.ps b2plot.png

The -density option increases the bitmap resolution. If you encounter the error "not authorized", find the file /etc/ImageMagick-6/policy.xml and change the line

<policy domain="coder" rights="none" pattern="PS" />

to

<policy domain="coder" rights="read|write" pattern="PS" />

Simple 1D plots can be converted to values using the write command, for example:

echo "phys sx xpar write 19 f.x" | b2plot

The numbers can be found in b2pl.exe.dir/b2plot.write.

Other built-in post-processing routines

Checking the perpendicular transport coefficients

In SOLPS-ITER working environment and in the runfolder, run:

2d_profiles
xyplot < ke3da.last10 #electron heat conductivity
xyplot < ki3da.last10 #ion heat conductivity
xyplot < dn3da.last10 #particle diffusion

Alternatively, you can load the files with Python or something similar. They are human-readable and contain only the radial positions relative to the OMP separatrix and the diffusion coefficient values.

Alternatively, go straight to the source of the 2d_profiles command and load the b2time.nc file. Sample Python script:

import xarray as xr
b2time = xr.open_dataset('run_path/b2time.nc')
Dn = b2time['dn3da'].data[-1]
pl.plot(Dn)   # use the radial cell number as the x axis

As a note, forcing a pedestal in the temperature profile is primarily done through making an artificial transport barrier in the transport coefficients profile (a dip around the separatrix). In density, a "pedestal" forms more readily, so usually constant profiles can be used.

Converting between poloidal, parallel and perpendicular fluxes and flux densities

Plasma transport in the B2.5 code happens in two independent directions: poloidal (x) and radial (y). The plasma solution is symmetric in the toroidal coordinate (z). The underlying transport physics is, however, defined mostly in terms of transport parallel to the magnetic field (\(\parallel\)) and transport perpendicular to it (\(\perp\)). As a result, users of SOLPS must learn to convert between these coordinates.

For the sake of proper understanding, we will devote three boring paragraphs to nomenclature. Proper names of fluxes can be a mouthful, so they are often shortened and interchanged. Usually you can tell what is what by the context and units, but brevity comes at the cost is clarity (and sounding uneducated). If you know the difference between a flux and a flux density, energy flux and heat flux, and poloidal, parallel and perpendicular fluxes, you can skip the next three paragraphs.

A flux is "amount of something passing through a given boundary per second". There particle fluxes (unit particles/s or s^-1), momentum fluxes, or energy fluxes (symbol \(Q\), unit J/s or W). Anyone calling energy flux "power flux" does not understand the beauty of analogy. Power does not flow; energy does. A flux density is "amount of something passing through a unit area of boundary per second". There are particle flux densities (symbol \(\Gamma\), unit s^-1m^-2), momentum flux densities, or energy flux densities (symbol \(q\), unit W/m²). A flux is always tied to a specific boundary, as in "the energy flux impinging on the outer target is 97 kW". A flux density, on the other hand, is a standalone physical quantity, as in "peak divertor target parallel energy flux densities can reach 200 MW/m² in fusion reactors". When discussing physics, you're usually speaking in terms of flux densities, not fluxes.

Energy flux refers to any type of energy (thermal, potential, electric...) flowing somewhere. It is the proper term for total divertor energy flux, convective energy flux, potential ionisation energy flux, electric energy flux et cetera. Heat flux refers only to the transport of energy in the form of chaotic thermal movement (heat) along asymmetries in the velocity distribution function (conduction). Heat flux is synonymous for conductive energy flux. "Conductive heat flux" is superfluous. "Convective heat flux" is wrong. "Total divertor heat flux" is also wrong, though it can be referred to as the total divertor heat load, assuming that the energy impinging on the target is immediately converted into heat and this heat is what you're describing. I use the symbol \(q\) both for heat and energy flux densities, because I haven't found a better alternative.

Poloidal fluxes (e.g. poloidal energy flux \(q_{pol}\) or \(q_x\)) are more or less a construct of SOLPS (and other simulation codes). The real physical nature of this flux is transport along magnetic field lines, without the hindrance of Larmor rotation. However, as 2D codes will assume toroidal (z) symmetry, it's useful to have the other two coordinates perpendicular to the toroidal coordinate. "Poloidal physics" is just a projection of the parallel physics, with the advantages that the simulation domain is shorter and easier to grid, and the results are easier to visualise. Parallel fluxes (e.g. \(q_\parallel\)) flow along the magnetic field lines. They are the defining feature of SOL transport physics, as they provide a quick connection to a strong energy and particle sink (divertor target or limiter). Perpendicular fluxes (e.g. \(q_\perp\)) are, in the context of this section, parallel fluxes projected onto the normal of the divertor target (or limiter). (Disregard for a moment that perpendicular can also be used as a synonym for radial, as in "perpendicular anomalous diffusion coefficients" \(D_\perp\), and forget all about the binormal \(\perp\) direction.) Conceptually, perpendicular fluxes only make sense in direct relation to the target, unlike parallel fluxes which are defined throughout the plasma volume. Perpendicular energy flux densities are the relevant quantity to the survival of the plasma-facing components. If an engineer is asking you how to design the divertor cooling systems, they are asking for the total perpendicular energy flux density. Even though the poloidal, parallel and perpendicular fluxes are all projections of the same phenomenon (parallel transport), they are applied in different fields (respectively: transport codes, SOL physics, and plasma-wall interaction), which subtly distinguishes them from one another.

To wrap up the last three paragraphs: You will often hear "heat flux" instead or "parallel energy flux density". The reason is prosaic - it's shorter. It's up to you to understand from the context if it's flux or flux density, heat or just energy, and what its direction is. This is fine in everyday communication, but before we proceed with converting between them, we need to clarify what we are talking about.

The following describes the calculation of the parallel and perpendicular energy flux density from an arbitrary poloidal energy flux \(Q_x\) [W], for example the total poloidal energy flux carried by the plasma fhtx. The constituents of the total target heat load will be described elsewhere.

1. Parallel energy flux density \(q_\parallel\) [W/m²]:

   \(q_\parallel = \dfrac{Q_x}{b_x s_{x,perp}}\), where \(b_x\) and \(s_{x,perp}\) are defined below

   This is the energy flux density relevant to parallel SOL physics, such as heat flux limiting or the heat flux fall off length \(\lambda_q\) calculation.

2. Perpendicular energy flux density \(q_\perp\) [W/m²]:

   \(q_\perp = \dfrac{Q_x}{s_x} = q_\parallel . \sin \alpha\), where \(s_{x}\) and \(\alpha\) are defined below

   This is the energy flux density relevant to target heat loads.

The quantities in the above relations are:

1. Magnetic field pitch \(b_x = B_x/B\)

   This is the ratio between the poloidal magnetic field \(B_x\) [T] and the total magnetic field \(B\) [T]. It is also the sine of the angle between the parallel direction and the poloidal direction (the magnetic pitch angle, which is not the impact angle denoted as \(\alpha\) in this section).

   SOLPS-postproc: Available in sim.geometry['bb']. \(B_x\) is bb[0,:,:] and \(B\) is bb[3,:,:].

   b2plot: \(B_x\) is called bx and \(B\) is called bb.

2. "Poloidal area perpendicular to the x-direction" \(s_{x,perp}\) [m²]

   This is the geometric factor omnipresent in SOLPS-ITER equations: \(s_{x,perp} = \sqrt{g}/h_x\), where \(\sqrt{g}\) [m³] is the cell volume and \(h_x\) [m] is the cell length in the poloidal direction.

   SOLPS-postproc: Available in sim.geometry. The cell volume is sim.geometry['vol'] and the poloidal cell length is sim.geometry['hx'].

   b2plot: \(s_{x,perp}\) is called sxperp.

3. "Poloidal area of cell contact" \(s_x\) [m²]:

   This is the area of the target surface which is covered by the cell/flux tube.

   SOLPS-postproc: Available as sim.geometry['gs'][0,:,:].

   b2plot: \(s_x\) is called sx.

4. Divertor target impact angle \(\alpha\) [radians]

   This is the angle between the magnetic field line and the target surface. In other words, \(\alpha= \frac{\pi}{2}\) - (angle between the field line and the surface normal).

   PLEQUE (equilibrium reconstruction handler): Can be calculated by evaluating the impact angle sine at the point of impact, \(\alpha\) = np.arcsin(eq.first_wall.impact_angle_sin()). (This will give you impact angles along the entire first wall.)

   SOLPS-postproc: Can be calculated using the previous three quantities stored in sim.geometry as \(\alpha = \mathrm{arcsin}(b_xs_{x,perp}/s_x)\).

   b2plot: Can be calculated as \(\alpha\) = arcsin(bx sxperp / sx bb). (Recall that \(b_x\) = bx/bb.)

If you'd like to flexibly convert between poloidal energy fluxes [W] and parallel/perpendicular energy flux densities [W/m²], we recommend either learning SOLPS-postproc, or plotting sx, sxperp bx/bb with b2plot and saving the results as a text file using the write command (see Questions and answers: Can I access b2plot data as numbers?).

A faster, out-of-the-box approach is to use the wlld post-processing command. Upon long and careful examination by Aleš Podolník, it has been concluded that not only does wlld correctly consider all the components of the energy flux impinging on the divertor targets, it also correctly converts them to the perpendicular energy flux density as \(q_\perp = Q_x/s_x\). If you don't want to use SOLPS-postproc, follow these steps:

1. Procure a converged, believable simulation and, in the SOLPS-ITER work environment inside the run directory, run:

   echo "1100 wlld" | b2plot
   

   This produces (among others) the files b2pl.exe.dir/ld_tg_i.dat and b2pl.exe.dir/ld_tg_o.dat, which contain perpendicular energy flux densities [W/m²] to the inner and outer divertor target, respectively. Remove them from b2pl.exe.dir/ and store them somewhere safe, as the next call of b2plot will overwrite the directory.

2. Export the inner target impact angle sines with b2plot:

   echo "phys bx sxperp m* sx m/ bb m/ write 0 f.y" | b2plot
   

   The part bx sxperp m* sx m/ bb m/ performs the calculation bx sxperp / sx bb. The write command saves the results into the text file b2pl.exe.dir/b2plot.write. The command 0 f.y picks only the radial profile at the inner target's first real plasma cell (the inner target ghost cell has the index -1). Remove the text file from b2pl.exe.dir/ and store it somewhere safe.

3. Export the outer target impact angle sines with b2plot:

   echo "phys bx sxperp m* sx m/ bb m/ write 83 f.y" | b2plot
   

   This assumes that you have followed the DivGeo tutorial and have 86 cells in the poloidal x direction, including the two ghost cells. If your nx is different, modify the 83 f.y command accordingly. Remove the text file from b2pl.exe.dir/ and store it somewhere safe.

4. Load the text files into your post-processing tool of choice and convert \(q_\perp\) into \(q_\parallel\) using the impact angle sine as \(q_\parallel = q_\perp / \sin \alpha\).

EIRENE neutral spectrum

Highly specific neutral-related data can be extracted from EIRENE during SOLPS-ITER or standalone runs using so-called tallies. Tallies are simply quantities calculated by adding up individual Monte Carlo results of an EIRENE run. Their definitions are generally given in the input.dat file and their values are typically printed in the standard output of EIRENE, which is typically redirected into run.log. In this example, we will show how to extract the neutral spectrum impinging on several vacuum vessel segments. Much information can be found in the EIRENE manualopen_in_new.

1. Define the surfaces where you want the tallies to be computed. These can be real vessel wall segments, completely virtual surfaces (such as the separatrix surface) or anything inbetween. The SOLPS-ITER input.dat file contains, by default, a complete list of vaccum vessel segments in block *** 3b. Data for additional surfaces. You can use these or define your own additional surfaces. A surface block looks like this:

   *   1 :   1
    2.00000E+00 1.00000E+00 1.00000E+00 1.00000E-05
        1     1     0     0     0     1     0     0     0    -1
   6.65934E+01 2.44670E+01-1.00000E+20 6.66467E+01 2.70939E+01 1.00000E+20
   

   Use the EIRENE manual and reverse engineering to figure out what the numbers mean. At the end, don't forget to edit the total number of additional surfaces, just under *** 3b. Data for additional surfaces.

2. Define how your tallies are calculated. In the example of neutral spectra, you need to edit block *** 10. Data for additional tallies and block ** 10f. Data for energy spectra. Refer to section 2.10.6 of the EIRENE manual.

   *** 10. Data for additional tallies
   0     0     0     0     0     8
   

   Here, 8 is number of spectrums that will be produced at the 8 (additional) surfaces which will be written out in block 10f. For one such surface, the corresponding lines are:

   ** 10f. Data for energy spectra
   11     1     1     1  -100     0     0
    2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00
   

   Here, 11 is the number of the surface where the tally is calculated, -100 is the number of bins in the spectrum, and 2.00000E+00 2.00000E+03 is the region of the considered energies. Write such two lines for each surface where you want the tally to be computed.

3. Tell EIRENE to actually print the tallies in block *** 11. Data for numerical and graphical output.

   *** 11. Data for numerical and graphical output
   FFFFF FFFFF FFFFF FFFFF
   TFFFF FFFFF FFFFF
        0     1
        1
       15     0     0     0     0    75
   TTFTT FFTFF FTFFF F
        1    37     1     1    87     1     0     0     0     0
   F                          0
   F                          0
   F                          0
   F                          0
   F                          0
   F                          0
   F                          0
   F                          0
    6.58893E+01 6.58893E+01 9.14404E+01-4.75189E+00 0.00000E+00
    1.00000E+01 1.00000E+01 0.00000E+00 0.00000E+00 1.00000E+01 1.00000E+01
    1.00000E+01 1.00000E+01 0.00000E+00
        0     0     0     0     0     0     0     0     0     0     2
        0
   

   The important number here is the 1 in line 0 1. This is the NSPCPR flag for printout of surface- or cell-based spectra, defined in block 10f. If NSPCPR > 0, the spectra are printed.

4. Run the code (coupled SOLPS-ITER or standalone EIRENE) and find the tally values in the EIRENE printout. This can be the standard shell output (what appears in the SOLPS work environment command line when you run b2run b2mn), run.log or run.log.last10 depending on the way you ran the code. Look for the "spectrum" keyword.

   SPECTRUM CALCULATED FOR ADDITIONAL SURFACE     15
   

   Visualise the spectrum using your favourite scientific software.

Change the number of EIRENE particles. To get sufficient statistical data (or very quick results), you can define the number of EIRENE particles to be launched in block 7 (should be checked with the EIRENE manual). Also, in block 1,

2     1    20   100     1     1     0     1     0

the number 100 is the number of particles to be launched.

Producing tallies with standalone EIRENE. Sven Wiesen recommends to run only EIRENE on frozen plasma background. This can be done by the command eirobj. 1. Create a new folder on the same level as your original SOLPS-ITER run. 2. Copy there the files fort.30, fort.31, fort.32, fort.33, fort.34, fort.35, fort.10, fort.11, fort.13, fort.15 and input.dat from your run folder. 3. Inside the case directory (on the same level as baserun), run correct_baserun_timestamps. 4. Run eirobj.

We have produced our spectra by running coupled SOLPS-ITER for several iterations.

Tallies of neutral fluxes. If you only need neutral fluxes (without the spectral resolution), you only need to modify only block 11. Refer to the EIRENE manual.

Example of modified input.dat. This will produce spectra at 8 vacuum vessel segments.

Block 1: change the NFILE flag, controlling where data is read from and which output files will be saved.

    *** 1. Data for operating mode
    # previously:    -1     1    50   406     1     1     0     1     0
        -1     1    50   100     1     1     0     1     0

Block 10: change the NADVI and NADSPC flags, controlling the number of tallies of each kind to compute.

    *** 10. Data for additional tallies
    # previously: 8     0     0     0     0     0
    0     0     0     0     0     6

Block 10f: set how to calculate the spectra.

    ** 10f. Data for energy spectra
        15     1     1     1  -100     0     0
     2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00
        16     1     1     1  -100     0     0
     2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00
        17     1     1     1  -100     0     0
     2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00
        43     1     1     1  -100     0     0
     2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00
        58     1     1     1  -100     0     0
     2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00
        59     1     1     1  -100     0     0
     2.00000E+00 2.00000E+03 0.00000E+00 0.10000E+00 6.66000E+02 1.00000E+00

Block 11: turn on the spectrum printout. Previously, the first half of block 11 read:

    *** 11. Data for numerical and graphical output
    FFFFF FFFFF FFFFF FFFFF FFFFF FFFFF FFFFF
    TFFFF FFF
         0
         0
    FFFFF FFFFF FFFFF F

Now, the first half of block 11 reads:

    *** 11. Data for numerical and graphical output
    FFFFF FFFFF FFFFF FFFFF
    TFFFF FFFFF FFFFF
         0     1
         6
        15     0     0     0     0    75
        16     0     0     0     0    75
        17     0     0     0     0    75
        43     0     0     0     0    75
        58     0     0     0     0    75
        59     0     0     0     0    75
    TTFTT FFTFF FTFFF F

The second half of block 11 remains unchanged:

         1    37     1     1    87     1     0     0     0     0
    F                          0
    F                          0
    F                          0
    F                          0
    F                          0
    F                          0
    F                          0
    F                          0
     6.58893E+01 6.58893E+01 9.08073E+01-4.75189E+00 0.00000E+00
     1.00000E+01 1.00000E+01 0.00000E+00 0.00000E+00 1.00000E+01 1.00000E+01
     1.00000E+01 1.00000E+01 0.00000E+00
         0     0     0     0     0     0     0     0     0     0     2
         0

Matlab

Quoting email from Irina Borodkina, 9 April 2020:

Dear All,

I would like to describe a bit the output files of SOLPS-ITER run, uploaded in the directory /run.4MW. Maybe this information will be helpful.

The B2.5 grid (for charged species) is saved in free format in b2fgmtry, while the EIRENE grid is a bit more complicated but saved in free format in files fort.33-34-35. The plasma variables are instead saved in free format in b2fstate and b2fplasmf. The neutral quantities on the triangular grid are stored in the fort.46. All these files (b2fgmtry, b2fstate, b2fplasmf, fort.33-34-35) can be read by some nice Matlab routines that you can find in $SOLPSTOP/scripts/MatlabPostProcessing/IO, for example

read_triangle_mesh.m - This will load fort.33,fort.34, and fort.35

read_ft46.m - This will load fort.46

read_b2fgmtry.m - This will load b2fgmtry

read_b2fstate.m - This will load b2fstate

Those routines are easily accessible and in there you can see which variables are being read. If some of the variable names are not clear you can look into the manual (SOLPS-ITER) in the chapter regarding b2plot or b2cdcv.

I also would like to ask to consider the uploaded files like preliminary results for technical testing. When you will be ready to import them to cherab, please let me know, and I will provide you with the best results I had by that time.

I wish a Happy Easter to all!

Cheers,

Irina

Python - overview

Currently (June 2024), I know of seven code packages which provide access to SOLPS output through Python:

• The SOLPSpy package, developed at ASDEX-Upgrade by Ivan Paradela Perez and Ferdinand Hitzler. (Don't worry if you can't open the link - further in the text we go through the process of how to get there.) Its installation and usage are described below.
• The Quixote package, a new, better and publicly available version of SOLPSpy developed by Ivan Paradela Perez. I've never tried using Quixote yet, but it has nice public documentation (unlike SOLPSpy).
• The Aurora package, which I have not tried personally but seems to have the same functionality as Cherab.
• The Cherab package, primarily written for tokamak spectrometry, can also load and process SOLPS-ITER output. Cherab is accessible at IPP Prague servers as a module; more information can be given by Matěj Tomeš (tomes@ipp.cas.cz).
• The solps-postproc package, built on top of Cherab and developed at IPP Prague by Jakub Seidl and Matěj Tomeš. This package is still under construction, but it promises not only accessing SOLPS output but constructing a "physical reality" from the discretised code results using interpolation and nomenclature independent of SOLPS particulars.
• The ERO2.0 code, [Romazanov 2020] has built-in capabilities to load SOLPS-ITER data.
• The EireneX package can load and post-process EIRENE data. As the package GitHub is currently unavailable, it is unknown whether it still maintained.

SOLPSpy

Until succeeded by Quixote, SOLPSpy was developed at ASDEX-Upgrade. Its source code is stored at AUG GitLab. This means that people not working directly on the AUG work stations can't access its home page. To see the page, use the VPN connection described in Remote access.

The master branch has not been updated in ever, so use the develop branch instead (dev). There are two ways to get the source code: clone it using git or download it.

Clone the SOLPSpy source code using a terminal

This option is available if you want to download SOLPSpy onto your AUG workstation. First, log into your AUG workstation via remote desktop or via SSH (see the same section of Remote access), go into the folder where you want to install SOLPSpy and clone the git repository.

cd ~/Documents
git clone ssh://git@gitlab.aug.ipp.mpg.de:2222/ipar/solpspy.git

After that's done, check the available branches:

cd solpspy
git branch -a

The develop branch (called dev) should be among the output. This is the branch you want to follow, as the master branch is not kept up very well. To switch there:

git checkout dev

To run SOLPSpy, you'll need to switch from the Solaris environment (saug23) to one of the AUG TOK Linux clusters (e.g. tok01). To do that:

rlogin tok01

Your username is filled in automatically (it's the same as your username on AUG intranet), your password is the same as the one for AUG intranet. If you are denied access after writing this password, ask David Coster to ask the appropriate person to give you access rights to the TOK cluster. Assuming you've had success, you'll need to load some modules now.

module avail

This will list all available modules. The ones we want are MDSplus (SOLPS-ITER output data format) and Anaconda with Python 2. To load the two modules:

module load mdsplus
module load anaconda/2/2019.03

Finally, tell Python where to look for the SOLPSpy module:

setenv PYTHONPATH /afs/ipp/home/k/katej/Documents

Download the SOLPSpy source code manually

If you aren't working on an AUG workstation, this is the way to go. Connect to IPP Garching VPN (the tutorial is in the Remote access document) and go to the SOLPSpy homepage. (If you don't have access to the AUG network, ask one of the SOLPSpy authors to send you the SOLPSpy source code instead.) On the homepage, switch to the dev branch (left red frame) and download the code (right red frame):

Extract the files from the archive and add the folder to your $PYTHONPATH:

export PYTHONPATH="${PYTHONPATH}:/path/to/solpspy"

To avoid having to change the $PYTHONPATH variable every time you want to run SOLPSpy, add the above line to your ~/.bashrc (or equivalent) file.

Meeting SOLPSpy dependencies

SOLPSpy is written in Python 2, so you need to load the according modules. Both on AUG workstations (after using rlogin) and on Soroban, Python (and iPython) is available via a module. To list all modules:

module avail

If the output of this command is too long to read, first find which shell type you are using:

echo $SHELL

If you're using bash, you can list modules whose name contains the word python in the following way:

module avail -t 2>&1 | grep -i python

If you're using tcsh, the same is done this way:

(module avail -t) |& grep -i python

Alternatively, look for the keyword anaconda and not python. In either case, choose carefully between the available modules and then load the one you need:

module load python/2.7-anaconda

To my knowledge SOLPSpy needs only the Python module; an MDSplus module is optional. This probably changes if you work with cases saved in the MDS+ format.

Be careful not to load two versions of Python at once. To list all loaded modules:

module list

To unload an unwanted module:

module unload python/3.7-anaconda

Using SOLPSpy

Once you have a folder with the SOLPSpy source code, no further installation is required. Just add the folder to your $PYTHONPATH (explained above), run iPython

ipython

and import SOLPSpy:

import solpspy

If it works, great! If it doesn't, that makes me sad.

Now we'll try loading SOLPS-ITER output. Get the path your run/ directory and call:

rundir = solpspy.RundirData('/path/to/the/run/folder')

If your run is not located on the same server, follow the instructions above.

Your new variable rundir contains all the data from your b2fstate file and more as attributes. They are all listed in this table. A few major ones have docstrings, the rest can be found in the manual or in the code source. For instance, jxi and jxa are midplane coordinates.

Since SOLPSpy runs inside a Python 2 terminal and my CDB data processing routines are in Python 3, I need to pass the rundir data from one terminal to another somehow. At the moment I'm doing this by constructing a Pandas dataframe from selected data, saving it as a pickle and loading the file in a Python 3 terminal.

Cherab

To use Cherab on COMPASS, follow these steps.

1. Open a terminal where your SOLPS-ITER simulations are stored (here Soroban front, which has links to node 06, my private SOLPS-ITER installation).

   ssh -X username@soroban   #add @tok.ipp.cas.cz if needed
   

2. Import the cherab/dev module.

   module load cherab/dev
   

   If this fails, purge all modules using module purge and try again. In the long run, configure your ~/.bashrc file, which contains the instructions performed every time you open a terminal, to deal with the module conflicts.

3. Open an iPython terminal.

   ipython
   

4. Import the function for simulation loading and load your SOLPS-ITER simulation of choice using its path.

   from cherab.solps import load_solps_from_raw_output
   sim = load_solps_from_raw_output(path_to_run_directory)
   

   If Cherab complains about missing file b2fgmtry or b2fplasmf, follow instructions below.

Browse the sim variable properties and methods. It gives you access to selected simulation results (\(T_e\), \(T_i\), \(u_a\), \(n_a\) etc.) through 2D NumPy arrays, which you can use for profile plotting. There are also beautifully quick and responsive 2D plots (much better alternative to b2plot).

ax = pl.gca()
ax.set_title('Electron temperature [eV]')
sim.mesh.plot_quadrangle_mesh(solps_data=sim.electron_temperature, ax=ax)
ax.get_figure().colorbar(ax.collections[0], aspect=40, label='$T_e$ [eV]')

Providing up-to-date files to Cherab

For full functionality, Cherab demands two files which may not be present in the run folder: b2fgmtry and b2fplasmf.

The geometry file b2fgmtry is most likely located in baserun if you ran b2run b2ag there. You can make a soft link to it from the run directory using:

ln -s ../baserun/b2fgmtry b2fgmtry

The plasma solution file b2fplasmf is a human-readable version of b2fplasma, which contains all the information of b2fstate and additional data calculated from it. If you followed the instruction in the manual, you would simply create it in the run directory using:

b2run b2uf

This command has backfired (deleted the plasma solution) so many times that it inspired Honza to write the tutorial Common pitfalls: The danger of timestamps. Refer to it.

Warning: b2fplasmf is not updated automatically. This means you must run b2run b2uf each time you run your simulation, or each time before you use Cherab. To save you work, here's a script by Honza:

#!/bin/bash

if [ -z "$SOLPSTOP" ]; then
  echo >&2 "Please init SOLPS-ITER environment first (source setup.csh)"
  exit 1
fi

correct_b2yt_timestamps
touch b2mn.prt b2fparam b2fstate b2fmovie b2ftrace b2ftrack b2fplasma

if (b2run -n b2uf | grep -qE 'b2mn.exe'); then
  echo >&2 "Simulation seems to be unfinished or crashed"
  exit 1
fi

b2run b2uf

Save this as a text file somewhere in your SOLPS installation (the case folder has done well for me). Mark the file as executable:

chmod +x create_b2fplasmf

Then, in your run directory, you can run the script using:

../create_b2fplasmf

...and wait a few minutes.

SOLPS-postproc

Thanks to Honza Hečko and Jakub Seidl, solps-postproc is available as a COMPASS module. This tutorial shows how to get started in using it. If you're interested in developing solps-postproc instead, refer to the following section.

How to start using SOLPS-postproc:

1. Log into the Soroban front and load the SOLPS module.

   ssh -X soroban
   module load python/3.9-solps
   

2. As prompted by the previous command, run:

   solps_to_jupyter
   

   This will prepare the SOLPS environment for Jupyter, so you can use SOLPS-postproc in your Jupyter notebooks. You only need to do this once.

3. If prompted to do so by the module load command, run:

   cherab_link_openadas
   

   This is a hot fix of OpenADAS atomic data, which is missing from this Cherab installation. You only need to do this once.

4. Browse the example Jupyter notebooks and download the first example notebook to your machine. Mark where you downloaded it, so you can find it in the next step.

5. Enter the IPP Jupyter Hub on the web address https://jupyterhub.tok.ipp.cas.cz. You can do this from your local machine, as long as you're in the IPP internal network (e.g. the wi-fi called COMPASS) or connected to the IPP Prague VPN.

6. Log in with your LDAP credentials (for logging in to IPP work stations, COMPASS logbook etc.).

7. For PBS queue, choose workstations.

8. Press Start at the bottom of the page.

9. Find and open the example notebook you downloaded earlier. Play around with it. Reverse engineer. Have fun. If something's missing or doesn't work, create an issue in the SOLPS-postproc GitLab page.

If solps-postproc ever complains about missing files b2fgmtry or b2fplasmf, refer to the above section Providing up-to-date files to Cherab. Also note that b2fplasmf must be manually updated prior to loading the simulation.

Becoming a SOLPS-postproc developer

The developers of solps-postproc, Jakub Seidl and Matěj Tomeš, take their work seriously and adhere to high standards. To take part in developing the package:

1. Set up an IPP GitLab account (some info is on the COMPASS wiki). The login credentials are the same as your IPP username and password.

2. Ask Jakub Seidl to add you to the DMS group as a Developer.

3. Log into Soroban front, load the newest Python environment and install Poetry.

   ssh -X username@soroban
   module load python/3.10-anaconda-2023.02
   curl -sSL https://install.python-poetry.org | python3 -
   

4. So your system knows where to look for Poetry, add ~/.local/bin to your $PATH:

   echo "export PATH=${PATH}:${HOME}/.local/bin" >> ~/.bashrc && source ~/.bashrc
   

5. Go to the solps-postproc directory and configure Poetry to create virtual environments in the project directory instead of tmpfs. (While the use of tmpfs could be arguably faster, it is local to each machine in the COMPASS network and often quite limited in size.)

   cd /my/favourite/folder/solps-postproc
   poetry config virtualenvs.in-project true
   

6. Disable keyring so Poetry can compile the virtual environment.

   export PYTHON_KEYRING_BACKEND=keyring.backends.null.Keyring
   

7. Build the virtual environment for developing solps-postproc with Poetry.

   module purge
   poetry install --with dev
   

8. Load the Python version which you will use to work with solps-postproc and tell Poetry to use it.

   module load python/3.9-solps
   poetry env use $(which python)
   

9. Check that you can access the SOLPS-postproc source code, then clone the solps-postproc source code into a director of your choice.

   cd /my/favourite/directory
   git clone git@repo.tok.ipp.cas.cz:dms/solps-postproc.git
   

   If you get an error message saying Please make sure you have the correct access rights and the repository exists., make sure you have copied your public SSH key to the IPP Prague GitLab. Follow this tutorial for instructions, especially sections See if you have an existing SSH key pair and Add an SSH key to your GitLab account.

10. Add the package location to your Pythonpath, so Python knows where to look for it.

    export PYTHONPATH=/my/favourite/folder/solps-postproc/:$PYTHONPATH
    

    Additionally, find the file ~/.bshrc in your home folder and paste the above line into it. This will execute the line every time you open a terminal from now on.

Developing SOLPS-postproc

Refer also to The SOLPS-ITER codebase: Git for dummies.

First, complete the instructions in the previous section. You only need to do this once.

To edit the code and manage Git, I recommend Visual Studio Code and GitAhead.

To contribute to the Git repository of solps-postproc:

1. Upgrade your local repository to the latest develop version. If needed, do a complete wipe (delete the directory and clone it anew).

2. In GitLab, Repository -> Branches -> Create new branch, based on the develop branch. Start its name with feature/, bugfix/, hotfix/, or test/ (usually feature/). More info on naming branches and commits is here.

3. Pull this branch to your local repository and checkout to it. If you're using GitAhead, check "Create detached HEAD", then when prompted, create a local branch of the same name.

4. Write your changes. Use the module load python/3.9-solps environment to run and interactively test solps-postproc.

5. Fix formatting.

   module purge
   cd /my/favourite/folder/solps-postproc
   scripts/format.sh
   

6. Diagnose syntactical and stylistic problems. Then correct them.

   scripts/lint.sh
   

7. Run tests. If errors are encountered, fix them and run tests again.

   scripts/test.sh
   

8. Group your changes into logical commits and commit them to your branch.

9. Push the branch to the remote repository.

10. In GitLab, create a Merge Request. Assign someone involved (probably Jakub Seidl) as the reviewer. If you are done in this line of work, request that the branch be deleted after merging.

11. Resolve any issues the reviewer might have. Keep committing to the same branch until the merge request is approved.

12. Once your branch has been merged into develop, pull the newest version and checkout to it.

Using solps-postproc in Jupyter notebooks

Jupyter notebooks allow writing code in a story-like form. To run solps-postproc in a Jupyter noteobok:

1. In the IPP network (or with active VPN), visit the IPP Jupyter Hub.

2. Log in with your LDAP credentials (for logging in to IPP work stations, COMPASS logbook etc.).

3. For PBS queue, choose workstations.

4. Press Start at the bottom of the page.

5. Find and open your notebook.

6. In the top bar, choose Kernel - Change kernel - SOLPS [Python 3.9].