Questions and answers

This document records the questions we've had about SOLPS-ITER and their answers. It is split into sections b2plot, Input files and boundary conditions, Post-processing SOLPS-ITER results, Running, restarting and transferring SOLPS-ITER runs, Convergence, Physics considerations in SOLPS-ITER and Miscellaneous.

Unanswered questions

The manual mentions a transfer_solps-iter_runs utility. However, using it at Marconi Gateway leaves me with sshpass: command not found and I don't think I can install it without superuser rights (which I don't have).

If realistic plasma potential is only reached in runs with drifts, is there any physical meaning to it in runs without drifts?

To SOLPSpy creators/users: Not all the rundir attributes are in the rundir class; where is te for example?

Does parallel distance from OMP in cells in PFR have any physical meaning? Plasma pressure goes all wonky at the end of the flux tube...

Note to self: Investigate echo "a4p 0. rmin vesl outl sq phys 0 6563 10 nspec same logf 3 ndec surf 'inside.chords' chpl" | b2plot &.

What is the correct value of the ion sheath heat transmission coefficient $\gamma_i$, 5/2 or 3/2?

Why is the kinetic energy flux to the target multiplied by zero in non-drift runs?

Formulas for some boundary conditions (like BCENE=12) are in the manual. What about the others? Only in-code?

Decoding wlld output: Why is Whtpl some 30 % higher than (fhe+fhi)/sx? Is there any other way to know than to delve into the wlld source code?

Is BCENE/BCENI=12 physically as high-quality as =3? It's much easier to interpret...

According to the manual, ns is "the number of atomic species in the calculation". Does this mean B2.5 neutrals don't include molecules?

Answered questions

`b2plot`

Should I take the time to learn `b2plot`?

According to Lisa Sytova, b2plot is not worth your time, with a few exceptions. One of them is the wlld command. For instance,

echo '1101 wlld' | b2plot

produces a folder with target heat fluxes, temperatures, densities and all that good stuff. It can also perform line-of-sight integration of individual species radiation (see Comparing experiment to model: Bolometry).

In my own experience, b2plot (and 2dt and resco and all that jazz) is handy while running a simulation. However, to produce any physics results, one's own scripts vastly outperform b2plot in most tasks. Even the mentioned line-of-sight integration can be done with software such as Cherab, and when you plot your own heat flux traces at least you understand what the contributions are.

Can I make "scripts" for `b2plot`?

Yes. Refer to Processing SOLPS-ITER output: b2plot.

The output of `b2plot` is an ugly tiny picture, can I get it in a higher resolution?

Every call of b2plot creates a file called b2plot.ps in your run folder. To open it, use:

gv b2plot.ps

The file is in the PostScript format, which is vector graphics and its resolution is as high as you like it. But be careful! Another call of b2plot overwrites the b2plot.ps file, so if you want to keep the picture, you have to rename it or move it somewhere else. Also note that the b2plot.ps picture will only be produced after you close the preview window generated by the b2plot command (for instance by clicking on it).

How do I plot the ion density of a given species using `b2plot`?

echo "phys na 1 zsel surf" | b2plot

The variable na contains ion densities of all species, 1 zsel selects the species number 1 (probably D+ depending on the simulation).

How do I set plot boundaries in `b2plot`?

Your experience with the zmin command and others like might be that they crop the plotted data, not manipulate the axes themselves. The zmin, zmax etc. commands are, indeed, the way to go. The reason b2plot does not respond terribly well to them is because it's trying to set the axis limits in a way that preserves the aspect ratio and simultaneously puts nice numbers on the axes. As a result, getting the view you want is a matter of careful optimization. If everything else fails (or you don't feel like waiting), turn the unity aspect ratio off by the aspc keyword and fit the aspect ratio by the eye.

The `b2plot` axes don't display properly. What do I do?

To fix this without closing the figure, change the size of the figure by dragging its frame. Gnuplot will start behaving again.

How do I label the axes with `b2plot`?

Upgrade to a SOLPS-ITER version released in 2021 and use the xlab and ylab commands.

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

When I run `b2plot`, it says `gnuplot: unable to parse 'BACKGROUND'. Using black` and the image background is black.

Try a combination of these:

module load gnuplot
module unload gnuplot

`2dt tesepm` produces a messy line.

Check the switch b2mndr_stim in your b2mn.dat file. If it's set to a negative value, the simulation results are appended to the previous runs. You might have restarted the run many times, resulting in jumps backwards in time.

How do I set colourbar extent not common for all plots?

When plotting several similar quantities together, b2plot automatically gives them all the same colour bar extent, rendering them unreadable. To set its extent individually for every plot, use:

echo "phys na 0.0 fmin 0.0 fmax surf" | b2plot

By setting the colour bar extent to 0-0, b2plot stops to think and uses a different extent for each plot.

Can I access `b2plot` data as numbers?

Yes! 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.

The `b2plot` picture is tiny and I can't see the details. Can I zoom in?

After you close the window b2plot creates, a file called b2plot.ps should appear in your run folder. This contains the vector version of the b2plot picture, which you can zoom in as you like.

Input files and boundary conditions

What is the default for boundary conditions? `b2ah.dat` ("standard" BCs) or `b2.xxx.parameters` ("physics" BCs)?

b2ah is the default because if the switches b2stbc_boundary_namelist, b2stbr_neutrals_namelist and b2tqna_transport_namelist are not specified in b2mn.dat, their value is assumed to be 0, meaning that b2ah.dat is used. However, this doesn't mean that you should use b2ah.dat for specifying boundary conditions.

How should the time step be chosen?

The time step is controlled by the line

'b2mndr_dtim' ' 1.00E-04'

in the b2mn.dat file. The bigger the time step, the bigger the changes that happen with each iteration. This can, on one hand, mean that the simulation converges faster, but if you choose too large a time step, the simulation can also get out of hand and crash. For a pure deuterium case you can start with 10^-4 s and then vary the time step in several different runs to check the result. Generally, the more complex your simulation (more impurities, drifts), the smaller time step you'll need. Refer also to Common pitfalls: Reduce the time step.

Where can I switch between EIRENE and standalone B2.5?

In b2mn.dat, find the lines:

'b2mndr_eirene' '1'
'b2mndr_rescale_neutrals_sources' '1.0e-10'

The first one switches EIRENE on and off (here it's on), the second one scales all neutral sources (see chapter A.4 of the manual). To turn EIRENE off, change these lines to:

'b2mndr_eirene' '0'
'b2mndr_rescale_neutrals_sources' '1.0'

Leave the b2mndr_rescale_neutrals switch as it is (at 1.0). Additionally, it is recommended to set the environment variable

setenv STAND_ALONE yes

Why are `b2mn.dat` and `b2.transport.parameters.stencil` not automatically generated?

It's on purpose. Xavier doesn't want people to blindly use automatically generated boundary conditions, so SOLPS-ITER only automatically generates b2ah.dat, b2.boundary.parameters.stencil and b2.neutrals.parameters.stencil. Any further input files (such as b2.numerics.parameters or b2.transport inputfile) have to be procured by the user, for instance by... copying them from the manual. (Sorry Xavier, but this was a lost battle.)

I copied the `b2.transport.parameters` file from the manual. What do all the 12's stand for?

The example file reads:

&transport
flag_dna=1, parm_dna=12*0.4,
flag_dpa=0, parm_dpa=12*0,
flag_vla=0, parm_vla=12*0,
flag_vsa=0, parm_vsa=12*0,
flag_hci=1, parm_hci=12*1.6,
flag_hce=1, parm_hce=1.6,
flag_sig=0, parm_sig=0,
flag_alf=0, parm_alf=0,
/

Here, 12* means "12 times the following value". It's handy when you have a lot of plasma species (you model a lot of impurities) and you don't want to clutter the file too much. This notation can be used in nearly all input files. As a note, in this particular simulation (D+C+He) there are 12 species including neutrals. When you use this file, make sure to adjust it according to the number of species in your simulation.

DivGeo and `b2.boundary.parameters` both define input power - which is the alpha male?

Actually, DivGeo defines the radiated power and b2.boundary.parameters define the ion and electron power crossing the separatrix. They're not the same.

Can I make comments in the SOLPS-ITER input files?

In some, yes. For instance anything that is written into b2mn.dat and does not follow the normal switch structure

'name_of_switch' 'its_value'

is ignored. So you can write comments there in any way you wish (starting the line with #, %, \<!-- …), the code will just ignore the extra text.

What's the sheath heat transmission coefficient value and how can I find it out?

Refer to eq. (14) in [Kotov & Reiter, PPCF 2009] and the B2.5 boundary conditions (given in the manual). For instance, if you're using BCENE and BCENI=3, you're specifying the electron and ion heat transmission coefficient value. Their role in B2.5 equations doesn't exactly correspond to what Kotov & Reiter write, but there is some qualitative agreement (they're both factors to something that looks like convective heat flux, proportional to density, temperature and flow velocity). In the end, there are two ways to calculate the total heat transmission coefficient (the one they use in probe measurements):

Calculate it according to equation (14) in Kotov & Reiter.
Get the total heat flux (fhi+fhe)/sx from SOLPS-ITER output and divide it by $e n_e T_e c_s$.

So far comparison of the two methods seems to yield a factor-of-two agreement based on the run.

What is the difference between the boundary conditions `BCENE=3` and `BCENE=12`?

In the manual, BCENE=3 reads "sheath conditions, electron energy transmission, ENEPAR(,1) specifies an additional contribution to the energy transmission coefficient in addition to that of the potential difference [the sound speed used depends on settings of MOMPAR(,ISMAIN,2)]". Conversely, BCENE=12 reads "sheath conditions, electron energy transmission coefficient, ENEPAR(,1) specifies an energy transmission factor, $\delta_e$ in $Q_e = \delta_e\Gamma_{e}T_{e}$". The difference is that while using BCENE and BCENI=12, Kotov & Reiter's formula for the sheath transmission coefficient actually works.

How do I give the diffusion coefficient a radial profile?

The terse and inevitable answer is, read the manual, section 3.5.5 b2.transport.inputfile and go through everything it refers to. The easy but incomplete answer is:

Create the text file b2.transport.inputfile in your run folder.
Example file content (you'll have to read the manual, sorry):
```
&TRANSPORT
ndata(1, 1 , 1 )= 3 ,
tdata(1, 1 , 1 , 1 )= -4.00E-02 , tdata (2, 1, 1, 1)= 1.00 ,
tdata(1, 2 , 1 , 1 )=  0.00E-02 , tdata (2, 2, 1, 1)= 1.00 ,
tdata(1, 3 , 1 , 1 )=  4.00E-02 , tdata (2, 3, 1, 1)= 4.00 ,
no_pflux=.true.
/
```
This translates to: "Create the diffusion coefficient $D_\perp^{na}$ profile in ion species with index 1 by linearly interpolating the following three points. 4 cm radially inside the OMP separatrix: 1 m$^2$s$^{-1}$. At the separatrix: 1 m$^2$s$^{-1}$. 4 cm outside the separatrix: 4 m$^2$s$^{-1}$."
Add the switch to your b2mn.dat:
```
'b2tqna_inputfile'          '1'
```
As a result, b2.transport.inputfile will override whatever is specified in b2.transport.parameters, b2mn.dat and b2ah.dat.
Perform a test run. Save your current b2fstate aside, set the simulation time to 60 s in b2mn.dat and see if the simulation crashes complaining about b2.transport.inputfile in run.log. On the first try, I usually discover I made a mistake in the indices.
After running your simulation for a bit, check the diffusion coefficient profiles using:
```
2d_profiles
xyplot < ke3da.last10 #electron heat conductivity
xyplot < ki3da.last10 #ion heat conductivity
xyplot < dn3da.last10 #particle diffusion
```
If the profiles don't seem right, double-check the indices of ndata and tdata.

What happens to B2.5 ions which reach the North or PFR South boundary?

(Many thanks to Xavier Bonnin, who answered this question on the SOLPS Slack, June 2024, channel #beginner-advice.) Short answer: they are neutralised and reflected back into the B2.5 computation domain, where their journey is followed by EIRENE. Long answer:

How EIRENE treats the N/S boundaries: In SOLPS-ITER (and SOLPS5.1), all particles (atoms, molecules, and molecular ions) followed by Eirene see the N/S boundaries from the B2.5 grid as transparent "switching" surfaces, but nothing else. They pass through them, in either direction, undisturbed and with their velocity vectors unchanged. They continue to travel until they reach the actual walls that are defined in block 3B the Eirene input file.
How B2.5 treats the N/S boundaries: The plasma ions followed by B2.5 are passed to Eirene and made to neutralize and reflect (and even possibly cause sputtering) on that virtual boundary as if it were a solid wall. The details follow the specifications from the surface block given for that boundary in block 3A of the Eirene input file for the corresponding non-standard surface. These ions are assumed to hit the virtual wall with an energy of $2 T_i$, with no sheath acceleration. Depending on the surface specification, a fraction RECYCF will be fast-reflected from the wall at this $2 T_i$ energy, times the energy reflection coefficient of that surface from the TRIM model. The fraction RECYCT will be thermally recycled at the temperature of the wall, as molecules in case of hydrogen isotopes.
Conclusion: The only particles that are allowed to cross the N/S boundary are the ones followed by Eirene. So a typical history will be: some ion flow is determined by B2.5 to cross the N/S boundary. This ion flow is reflected at the location of the boundary as if a wall was there. That reflection is handled by Eirene, which sends back a mixture of fast and slow neutrals into the domain covered by both the B2.5 and Eirene grids. The neutrals are then followed by Eirene. Some will be ionized back into the plasma, but others will endure collisions that will send them back roughly in the direction they came from, or have trajectories that make them cross the N/S boundary at another location. These neutrals will cross the N/S boundary unimpeded and continued to be followed by Eirene in the "void" regions between the B2.5 plasma grid and the actual vessel wall contour. That region is part of the domain covered by the Eirene triangular grid, so the trajectories can be followed accurately. Eventually these neutrals in the "void" regions will re-enter the plasma and get ionized, or find the pump get absorbed, or stick to the wall in places where the recycling coefficient is less than one.

Post-processing SOLPS-ITER results

How can I access SOLPS-ITER output?

Don't reinvent the wheel - use existing software and read Processing SOLPS-ITER output.

To get quantities at the target, should I use guard cell data (like Stangeby says) or last plasma cell data (like `wlld` does)?

Depends on which quantity. For fluxes, take the flux into the guard cell, for temperature and density take the last plasma cell. Also note that the fluxes are usually defined as the flux going into the cell on the left (x index smaller by 1), so the appropriate target flux indices will be different for the inner and outer target.

Are the SOLPSpy calculations the same as `b2plot`/SOLPS-ITER calculations? David Coster says one should not double the calculation...

Ivan tries his best to mimic the b2plot calculations, but he might be wrong in some cases.

`b2run b2uf` results in the error `error while loading shared libraries: libmkl_intel_lp64.so`.

This is very Soroban-specific, but go to $SOLPSTOP/SETUP/setup.csh.SOROBAN.ifort64 and add the line:

module load intel_oneapi

Paste this line into your SOLPS-ITER work environment and run b2run b2uf again.

I added plot zones to an existing run, but they don't show up in the output of `echo "1111 wlld" | b2plot`.

(Asked by Eric Emdee at the SOLPS Slack, and answered by Xavier Bonnin.) After declaring the plot zone in DivGeo, you have to Output the DivGeo file and re-run Uinp. But before you do re-run Uinp, you need to delete the existing files it had already created in your previous attempt. The Uinp screen output tells you that there are some files already found that it does not try to write again. The plot zones are, I believe, written out in b2.neutrals.parameters or b2.user.parameters (can't remember right now).

Running, restarting and transferring SOLPS-ITER runs

How do I check the status of my running SOLPS-ITER simulation?

Refer to Common pitfalls: How to recognise a simulation crash.

When I use the switch `b2mndr_elapsed` or `b2mndr_cpu` to set how long SOLPS-ITER should run, the code ends after 12 iterations (something like 15 seconds) regardless of the input value.

Set a value both for b2mndr_ntim and b2mndr_elapsed. Pick some ridiculously large value for the former (like 1,000,000) and your actual desired value (like 57600 seconds = 16 hours) for the latter. SOLPS-ITER will stop running when one of these conditions is violated.

Can cases from a different geometry be transferred as the initial state for another case?

Yes, using the command b2yt. Refer to the manual, section 3.14.

When I run `b2run` directly, I can see the files in the `run` folder updated as the simulation proceeds and I can track the progress using the `cpu` command. But when I use `itmsubmit`, the files remain unchanged and the `cpu` command thinks nothing is currently running.

The files are stored in the cluster memory for a while and only written to the hard disk when memory's low. Thus, it just takes a while for the simulation to show up in the run folder.

What are all the files which define a finished run? If I want to create a number of runs stemming for a single run, which files do I need to copy into the run folder?

The main hub of simulation results is b2fstate. This holds the basic quantities (densities, velocities, temperatures), from which SOLPS-ITER can calculate everything else. When you wish to seed a simulation with a previously converged solution with the same cell count and the (more or less) same ion species, all you need to do is:

cp origin_case/b2fstate target_case/b2fstati

(Substitute by scp if the target simulation is on a different server, refer to SOLPS-ITER user wisdom: Useful commands.)

For loading simulation results with post-processing software, you need a bit more files. Cherab and solps-postproc packages, for instance, require these files:

b2fstate
b2fplasmf  #up-to-date version!
b2fgmtry  # copied from baserun
input.dat
b2mn.dat

The list will vary depending on the software; some packages can also load the fort files, others may not need b2fplasmf etc. Refer to Processing SOLPS-ITER output for details.

Finally, the minimal variant of branching out a simulation is described in Remote access: Transferring SOLPS-ITER runs and Running SOLPS-ITER: Branching out a SOLPS-ITER run.

What is the initial solution of SOLPS-ITER simulations, so called "flat profiles"?

Excerpt from the manual, sec. 3.4:

The initial state being created by b2ai.dat will be a homogeneous isothermal static plasma, which is very far from a realistic plasma edge solution. It is therefore almost always better ans strongly recommended to start a case using an already converged case of the same topology and containing the same main species, even it does not have the same grid size, or corresponds to a different device, or if the species list is somewhat different. To convert such a solution to the parameters of the case at hand, the b2yt converter is provided.

There is one EIRENE runs for 14 B2.5 runs and it drives the residuals up every time. In what part of this cycle should one take the "final solution"? Just before the next EIRENE call, when the residuals are minimal? Or does it not matter?

The calculation always stops just before the EIRENE call. The periodic changes do, however, affect the plasma state solution uncertainties. To reduce the uncertainties, check the manual for EIRENE averaging schemes.

Does restarting the code affect the solution? For instance, is there a difference between running the code for 8 hours and running it for 4 hours and restarting it for another 4 hours?

There might be. In the words of Fabio Subba, paraphrased by David Tskhakaya, 2021:

Yes, the restarting of SOLPS-ITER is not perfect: there can be possible “problems” when restarting the run, or there can be some disagreement with nonstop runs; this is related to the fact that code do not (can not) save all the information on the physical and numerical state of the system (I assume). They show up rare, usually are explained and corresponding bugs are fixed, but there is no guarantee that there will be no new ones; so we have to leave with this.

Running `b2run b2yn` returns "no rule to make target". What do I do?

Create the file b2yn.dat in your run directory. An example is given in the manual, section 3.12 b2yn, display the progress of inner B2.5 iterations:

*idout    (int, char*80, free format)
0, 'cpte'

This file contains instructions what b2yn should actually do, so of course it won't start without it.

Convergence

See also Running SOLPS-ITER: Checking if the run has converged and Common pitfalls: Simulation divergence.

How can I tell if my converged simulation is "physical"?

Long story short, look at the values of temperature, density, heat fluxes etc. and compare them with experimental or expected values. Other than that, check your boundary conditions. If you had to use the anomalous transport coefficient D = 10⁶ m²s^-1, then something may be wrong. As a beginner, consult David Coster and let him have a look at the finished run. But there aren't really any other tricks besides these.

What is the average time for a simulation to converge?

People will tell you it's weeks to months, depending on time step, drifts, speed-up schemes, tokamak size etc. But so far my experience with pure D runs in COMPASS geometry, coupled B2.5+EIRENE, they take two hours to converge from "flat profiles" (initial state setup when you first create a run from scratch) and five minutes to converge from previous, similar simulations. Smaller machines are significantly faster to model, so COMPASS should be a breeze compared to ITER. Drift cases, for example, require a significantly smaller time step than non-drift cases (say, 10^-7 s instead of 10^-4 s), so once they are turned on the convergence time will skyrocket. As Sven has warned me on the MST topic 14 meeting in July 2020, the smaller the machine, the smaller the time step has to be in a drift case.

How do I speed up the convergence?

Refer to Common pitfalls: Make SOLPS-ITER run faster.

What is the meaning of residual plots produced by `resco/mo/po`? How are residuals calculated? What is their expected value? Why do my residuals seem to decrease exponentially for ~15 iterations and then spring up by two orders of magnitude?

Residuals come from the MHD equations that B2.5 solves. Generally, these equations have the form

(time derivative) + (divergence) = (sources)

or, after a simple rearrangement,

(time derivative) + (divergence) - (sources) = 0.

When you compute the right-hand side of the equation for each cell of the B2.5 grid in a converged run, you get something small, but generally non-zero. Residuals are what you get when you compute this small value for each cell, take the absolute value and sum them. Section C.10 of the manual contains the formula.

The value of residuals dramatically depends on what sort of case you're running. If it's a stand-alone B2.5 run (neutrals are modelled as fluids) in pure deuterium, residuals can get as low as 10^-7 (unit is unknown and, as explained later, inconsequential). If you're running a case with impurities and EIRENE neutrals, residuals may never climb below 30, no matter how many CPU hours you sacrifice. The thing is, the value of residuals themselves is inconsequential. What matters are two things:

Residuals and key plasma quantities don't change with further iterations. You can check this by running this command in the run/ folder:
```
gv b2time2d.ps
```
Here, gv is a programme that opens the file b2time2d.ps. If you see something like this:

...you know that further iterations are useless. Having arrived at a stable solution is a better indication of case convergence than any value of residuals.
The resulting plasma makes physical sense (verified by plotting, for instance, heat flux profiles on the target or the upstream temperature profile).

Finally, those weird spikes on your residuals? Those are EIRENE runs. They introduce a lot of Monte Carlo noise into the simulation, which is also why coupled runs (using B2.5 and EIRENE, that is, the entire SOLPS-ITER) tend to have larger residuals than stand-alone runs (using only B2.5). Check the values of

'b2mndt_nstg0' '1'
'b2mndt_nstg1' '1'
'b2mndt_nstg2' '15'

in the b2mn.dat file and compare them to how often you see a spike. (For description of what these switches mean, check section A.4 of the manual.)

Physics considerations in SOLPS-ITER

Why does changing the ion heat flux limiter change the electron heat flux?

It changes the ion temperature, which propagates to the electron temperature etc.

How important is exact separatrix position in SOLPS modelling?

Separatrix density is often used for feedback control in modelling, so knowing where the separatrix is important. On the other hand, separatrix temperature does not vary much across different regimes, so this value isn't of much importance. Rather than $T_{sep}$, precise knowledge of the power crossing the separatrix is needed. Finally, because there are larger density gradients in H-mode (due to pedestal), exact separatrix position is more crucial in H-mode than in L-mode.

How much should one trust the predictive capabilities of SOLPS-ITER (like in ITER or COMPASS-U cases)?

Efforts are underway to model future machines reliably on physics basis, but the codes tend to crash at the slightest disturbance. As for SOLPS-ITER, you have basically three options. You can make educated guesses as to the input parameters, for instance "I have these diffusion coefficients at ASDEX-U and I use the same ones for ITER". Another option are scaling studies - for instance, the SOL heat flux width in ITER is a hot topic. Finally, if you're really unsure which values to pick, pick all of them in turn and make a parameter scan. That way your engineers will at least know the order of magnitude.

Can you model 3D limiter geometry with SOLPS-ITER?

SOLPS-ITER is inherently 2D. That means 3D effects (a single lithium divertor tile, toroidally asymmetric limiters) cannot be modelled based on first principles. Depending on whether you suspect the 3D effects will play a large role or not, the solutions diverge. If the effect is strongly 3D (such as edge radiation modelling during RMPs), then switching to EM3-EIRENE may be more feasible. However, if the effect is only weak (such as limiter geometry in COMPASS-U when the plasma is limited and the B2 grid doesn't even come close to the limiters), then SOLPS-ITER may still be an adequate tool.

How is $P_{SOL}$ distributed along the separatrix?

The radial heat flux is defined by an anomalous heat diffusion coefficient (at the moment, constant poloidally and radially) and the temperature gradient. Due to the Shafranov shift, the temperature gradient is higher at LFS, so more power diffuses through there. To change this up, you can mess with the diffusion coefficient and, for instance, define ballooning power fluxes, but they don't make a great difference on the plasma state. They do, however, change the position of the stagnation point (where parallel plasma velocity is zero, toward neither target).

IR camera measures total heat flux (including radiation) to the target. What about the probes, do they only measure the charged particles contribution?

Yes, the probes only measure the charged particle power flux, that is, the electron and ion convective heat flux (conduction is assumed zero in the isothermal sheath) and the ion kinetic energy flux. The other parts of total target heat flux, recombination power, radiation power and incident neutrals power, can be negligible or significant depending wholly on the plasma conditions. Ask David Tskhakaya for more details or study the wlld output.

In the B2.5 equations, what is up with those tildes and "divergence-free" parts? Why does it influence convergence?

It's a long story. As [Rozhansky, 2001] describes in section 2.1, one of the terms in the particle flux used in the continuity equation is the diamagnetic flux, $$ V_\perp^{(dia)} = -\frac{1}{enB} \frac{1}{h_y} \frac{\partial (nT_i)}{\partial y} \hspace{1cm} \mathrm{and} \hspace{1cm} V_y^{(dia)} = \frac{B_z}{enB^2} \frac{1}{h_x} \frac{\partial (nT_i)}{\partial x}. $$ This term is stupid because (presumably) it contains derivatives of the ion pressure and this causes numerical instabilities when it's plugged inside the divergence. (That is, the divergence of the particle flux which forms one of the three terms in the continuity equation, next to the time derivative of the particle density and the particle sources and sinks.) Luckily, as [Bufferand, 2016] explains in formulas (11) and (12), the diamagnetic flux can be decomposed into two components: the so-called magnetization flux and the grad-B drift of the ion guiding centres. The magnetization flux divergence is, by definition, zero, and the second term in grad-B drift is usually negligible. This means that inside the divergence, the diamagnetic drift can be replaced by $$ \tilde{V}_\perp^{(dia)} = \frac{T_iB_z}{eb_z} \frac{1}{h_y} \frac{\partial}{\partial y} \left( \frac{1}{B^2} \right) \hspace{1cm} \mathrm{and} \hspace{1cm} \tilde{V}_y^{(dia)} = -\frac{T_iB_z}{e} \frac{1}{h_x} \frac{\partial}{\partial x} \left( \frac{1}{B^2} \right). $$ This is great news, because we no longer have to differentiate plasma parameters in every step. SOLPS-ITER is an electrostatic code, it treats the magnetic field as a given, so one can precalculate its derivative and lend the code a bit more stability. Or so they thought, until they tried modelling H-mode and all hell broke loose again. As a result, one more transformation was implemented, described in formulas (7-9) in [Rozhansky, 2009]. It uses the fact that it's the Pfirsch-Schlütter currents which close the circuit polarised by the grad-B drift. The grad-B drift is vertical and the PS currents flow along the field lines, so they don't have the same magnitude, but they do fulfil $\nabla( \textbf{j}_{grad-B} + \textbf{j}_{PS}) = 0$ due to charge conservation. As a result, the divergence of the grad-B drift is the same (except for the minus) as the divergence of the PS currents. And this is why you'll encounter terms with "P.Sch." written on them in the B2.5 equations.

As a note, the diamagnetic drift is, well, a drift. So all of this is only relevant to runs with drifts on. Also, I am still confused how to calculate the real particle fluxes. [Rozhansky, 2001] explicitly warns against using the grad-B drift to this end. The manual repeats this warning under equation (C.230). And then it goes ahead and uses the grad-B drift expression to calculate the real particle flux, neglecting the magnetization flux.

In power balance studies where should one take upstream, OMP or X-point?

Decide based on your simulation. Plot the relevant quantities (total heat flux and pressure) along the first SOL field line and check their maxima. In most cases you'll find them near the X-point, and thus you take the X-point as upstream.

How accurate is the parallel coordinate in SOLPS compared to PLEQUE?

Internally, SOLPS-ITER doesn't use the parallel coordinate at all. All calculations are made in the poloidal $x$ direction and projected using $b_x = B_x/B$. The only place where the parallel coordinate is used is in b2plot, switch xpar, which plots the parallel coordinate instead of the poloidal coordinate as the x axis. The parallel coordinate is computed in $SOLPSTOP/modules/B2.5/src/b2plot/fx.F, code block mentioning the keyword xparallel. This calculation is mirrored in SOLPSpy and SOLPS-postproc. The question is: since this calculation starts from the discretised B2.5 mesh, how accurate is it?

The figures show a B2.5 grid based on a COMPASS equilibrium, field lines traced using PLEQUE from the outer target cell centres (red dots) and the connection length from the OMP to the inner and outer target. Two observations can be made:

The poloidal cell rows follow PLEQUE field lines, but the field lines do not necessarily pass through the cell centres. This is most readily seen in the top left part of the equilibrium, but it is also true at the outer midplane (OMP). The takeaways are:
1. Poloidal cell rows in B2.5 truly constitute flux tubes, and the transport through them is purely parallel transport.
2. If, nevertheless, you want to trace field lines from the cell centres with PLEQUE, don't start from the OMP. You will follow a field line which does not pass through the poloidal cell row.
The discretised calculation of b2plot, SOLPSpy and SOLPS-postproc is, to a large degree, accurate. The only sizable deviation is seen near the separatrix, where the connection length goes to infinity near the X-point. This is a pretty good result for such coarse discretisation.

Where does the ion grad-B drift point?

Listing the ion grad-B drift direction is common courtesy while describing a SOLPS-ITER simulation. This is, in part, because the ion grad-B drift pointing toward the X-point facilitates access into H-mode. You can use PLEQUE to find its direction.

# Import the equilibrium-reading function from PLEQUE
from pleque.io.readers import read_geqdsk

# Load the equilibrium (must be in the EQDSK format, produced e.g. by EFIT)
eq = read_geqdsk('/path/to/equilibirum/eqdsk/file')

# Print the toroidal magnetic field at the magnetic axis
eq.B_tor(R=eq.magnetic_axis.R, Z=eq.magnetic_axis.Z)

Since the $\nabla B$ vector always points toward the main tokamak axis (passing through the donut hole), the direction of the ion grad-B drift $\textbf{B}\times \nabla B/B^2$ is determined purely by the direction of the magnetic field B. This is, in turn, the same as the direction of the toroidal magnetic field at the magnetic axis, calculated above. Since PLEQUE uses the COCOS field direction convention, negative $B_t$ means that the ion grad-B drift is pointing downward. To verify this, plot the field directions with eq.plot_geometry() and use the right-hand rule to calculate the direction of the vector multiplication.

Miscellaneous

How do I change the keyboard layout on AUG workstations?

Short answer, you don't change it, you get used to it. If you're accessing the IPP Garching Linux clusters remotely, consider using the ThinLinc Client instead of the Oracle Virtual Desktop Client, as this will not force the retarded keyboard layout on you.

The script `itmsubmit` fails, complaining something about wrong partition. I can't even print the help using `itmsubmit -h`.

Update SOLPS-ITER.

stop
solps-iter_update_develop

How do I change the password to my Gateway account?

Log into the account, open the command line and type:

passwd

The new password has to be at least 8 characters long.

How do I add a gas puff source inside the B2.5 grid region?

(Taken from the 36th SOLPS-ITER user forum debriefing.) The easiest method is to use DivGeo. Create a surface at the location of the gas puff inlet, just inside the grid, by creating two points and connecting them. Make sure that the surface normal for that surface points in the opposite direction to the direction in which you want the gas to be injected. Then, you must make that surface transparent to Eirene by setting, in "General Surface Data", the value for the "Surface Type" variable to "0". Additionally, it is best to set to 0 all parameters related to sputtering. Then, using the "Variables | Add | Gas puff" dialog box, specify the usual gas puff parameters (species, strength, etc…). After that, save your DG model and run Uinp to obtain the necessary Eirene and b2.neutrals.parameters.stencil input files.

The diagnostics are bollocks, what do I do with it? (Target $T_e$ especially…)

Check how you calculate the values. For instance, at the divertor targets correlations between $T_e$ and $n_e$ are significant, and so $q = \langle q \rangle = \langle T_e I_{sat} \rangle$ (raw data are multiplied, then averaged) gives different results than $\langle q \rangle = \langle T_e \rangle \langle I_{sat} \rangle$.
Look into the definitions of your quantities. For instance, a common assumption in calculating the sound speed is $T_e = T_i$. However, in most SOL conditions the ion temperature is higher, which can lead to an error. Make sure that your "total" quantities take all their components into account (electron + ion, static + dynamic, conductive + convective, parallel / poloidal etc.).
Look into your equilibrium. How did you produce the equilibrium file (standard EFIT reconstruction / Ondra Kovanda's better EFIT)? What input did you use (IPR coils / divertor Mirnov coils / flux loops)? If needed, build a new mesh and start over.

Where do I find the latest version of the manual?

At the ITER Sharepoint. Alternatively:

Download the latest version of SOLPS-ITER.

cd /my/favourite/directory
git clone ssh://git@git.iter.org/bnd/solps-iter.git

Go inside the cloned directory.
```
cd solps-iter
```
If you require a specific version of SOLPS-ITER, checkout to the branch.
```
git checkout carlis/feature/ad-wg
```

Update every module one more time.

git submodule init
git submodule update

Initiate the SOLPS work environment.
```
tcsh
source setup.csh
```
Compile the SOLPS-ITER manual.
```
gmake manual
```

You will find the result in $SOLPSTOP/docs/solps/solps.pdf.

When updating SOLPS-ITER, there are two options: `solps-iter_update(_develop)` and `git status; git fetch; git pull; git submodule update`. Is this it or are more steps needed?

The first option is complete while the second option isn't, as it doesn't compile the code. The git commands must be followed by make, make_all, make_all_mpi etc. For us beginners, therefore, the first option is the way to go. In summer 2020, the master branch was updated after 2 years of being woefully out-of-date, so the _develop isn't compulsory anymore.

How do I install a particular version of SOLPS-ITER?

On the second page of the manual, the code version from which it was built is written:

This manual is built from the following code repository versions:

SOLPS-ITER version : 3.0.7-41-g0c21b66

B2.5 version : 3.0.7-42-g5fe4159

Eirene version : 3.0.7-26-g9be0d5f

Carre version : 3.0.7-5-gef3eaef

DivGeo version : 3.0.7-6-g594e4d6

To make a SOLPS installation of this particular SOLPS version, make a new clone of the SOLPS installation, initialise and update the submodules:

git clone ssh://git@git.iter.org/bnd/solps-iter.git
cd solps-iter
git submodule init
git submodule update

Details will vary depending on where you're installing; see the installation manuals. Now's the time to switch to the desired version for SOLPS-ITER and its submodules.

git checkout 3.0.7-41-g0c21b66
cd modules/B2.5
git checkout 3.0.7-42-g5fe4159
cd ../Eirene
git checkout 3.0.7-26-g9be0d5f
cd ../Carre
git checkout 3.0.7-5-gef3eaef
cd ../DivGeo
git checkout 3.0.7-6-g594e4d6

If you've done this in the SOLPS-ITER work environment, exit it so the updated setup.csh can be executed. Finally, compile SOLPS according to the respective installation instructions.

exit  # only if you were already in the SOLPS work environment
tcsh
source setup.csh
gmake -j 4

What is a "stratum"?

A stratum (plural strata) is a neutral particle source of (in some regard) homogeneous properties. When sampling Monte Carlo neutrals whose trajectories it will follow, EIRENE splits its total neutral particle source into strata according to their origin - mostly "originating on a plasma boundary" or "volumetric neutral particle source". If it samples, say, 7000 neutral particles for each stratum, it achieves lower statistical variance than if it sampled from the total neutral particle source at random. This is because stratified sampling provides more representative samples of neutral particles in each EIRENE run.

More information can be found via Google (general understanding of stratified sampling in Monte Carlo techniques), in the EIRENE manual, sections 1.3.3.2 Stratified Source Sampling (particular implementation of stratified sampling in EIRENE) and 2.7 Input data for Initial Distribution of Test Particles (description of block 7 in input.dat, where strata description is given), and b2.neutrals.parameters along with the description of its switches.

For instance, the beginning of the b2.neutrals.parameters file in a D+C simulation reads:

 &NEUTRALS
 nstrai= 13,
 crcstra= 'W', 'W', 'E', 'E', 'S', 'S', 'S', 'S', 'N', 'N', 'V', 'V', 'T',
 rcpos= -1, -1, 84, 84, -1, -1, -1, -1, 36, 36,  0,  0,  0,
 rcstart=  0,  0,  0,  0,  0,  0, 66, 66,  0,  0,  0,  0,  0,
 rcend= 35, 35, 35, 35, 17, 17, 83, 83, 83, 83,  0,  0,  0,
 species_start=  0,  2,  0,  2,  0,  2,  0,  2,  0,  2,  0,  2,  0,
 species_end=  1,  8,  1,  8,  1,  8,  1,  8,  1,  8,  1,  8,  8,
 ...

It says there are 13 strata (nstrai): 10 associated with B2.5 grid boundaries (W, E, S and N), two volumetric sources (V) and one time-dependent source (T). The 10 boundary strata are, alternatively, deuterium (ion species 0-1) and carbon (ion species 2-8) neutral sources. The two volumetric particle sources also produce, respectively, neutral deuterium and carbon.