Feature blog aims to introduce additional important features of SOLPS-ITER which did not fit in the main step-by-step tutorials. Familiarity with the basic workflow is assumed.

Drifts terms and how to enable them

Quick checklist to enable drift terms (see below for details):

1. Enable drifts in b2mn.dat using the b2news_fac* switches.
2. Modify boundary conditions in b2.boundary.parameters to account for drifts.
3. Apply the speed-up schema in b2.numerics.parameters to improve numerical stability & computational speed.
4. Run and troubleshoot in case of crashes (dt, mesh, standalone, ...).
   • Decrease timestep, but no sense to go lower than 1e-7, 1e-8.
   • Adjust the fluid mesh (no large jumps in cell size, finer mesh, ...)
   • Start with simpler simulation (no Eirene, no impurities, no puffing, ...)

Motivation

By default, the plasma particle drift effects are not included in the B2.5 fluid equations and instead they are considered as one of the components of the mysterious "anomalous diffusion". The main reason is that including them is simply too computationally expensive and introduces numerical instabilities, possibly due to a use of different solver for the equations (unverified claim). Then again, under certain conditions (small machine, low magnetic field, ... see below), the drifts may play an important role and help you get far closer match to the experiment. This section explains how to enable them.

Pros	Cons
Better match with experiment (potentially)	Computationally expensive
Better accuracy of predictions (theoretically)	Emotionally expensive (numerical issues, results not guaranteed)
Dabbling with often ignored physics	Definitely not a magical solution for everything

Physics background - Drifts in a tokamak

construction This section is under construction.

1. Enable drifts with options in b2mn.dat

The drift terms are enabled using 3 groups of switches in b2mn.dat corresponding to 3 independent terms in the equations that can be activated:

• b2news_facdrift diamagnetic terms, ion inertia current
• b2news_facExB ExB drift terms
• b2news_facvis drift viscosity

Each group contains 4 switches that are used to control a multiplicative factor \(f\) in front of the term in the equation. The factor can be set to values in range \(f \in [0.0, 1.0]\), where \(0.0\) disables it and \(1.0\) enables it completely. The switches are named as the name of the group + one of the following suffixes:

• _start - initial value, \(f_0\)
• _target - target value, \(f_{target}\)
• _inc - multiplicative increment, \(f_{n+1} = A \cdot f_n\) on each timestep until \(f_n = f_{target}\)
• _dec - multiplicative decrement, applied when eqs. are not converging on that timestep

b2mn.dat

! A complete list, including the default values

'b2news_facdrift_start'  '0.0'
'b2news_facdrift_target' '0.0'
'b2news_facdrift_inc'    '1.0'
'b2news_facdrift_dec'    '0.0'

'b2news_facExB_start'    '0.0'
'b2news_facExB_target'   '0.0'
'b2news_facExB_inc'      '1.0'
'b2news_facExB_dec'      '0.0'

'b2news_facvis_start'    '0.0'
'b2news_facvis_target'   '0.0'
'b2news_facvis_inc'      '1.0'
'b2news_facvis_dec'      '0.0'

Possible approaches to setting the switches:

1. Constant profile \(f(t)=const.\)
2. set _inc = _dec = 1.0, then \(f_0 = f_n = f_{target}\) is kept during the whole run
3. increase _start = _target = 0.3 manually between runs (e.g. 0.3, 0.5, 0.7, 1.0)
4. Ramp-up profile, customize how steep
5. set _start = 0.1 (something small), _target = 1.0
6. set _inc = 1.5 (something larger than 1)
7. set _dec = 0.25 (something smaller than 1) as a fallback in case of divergence

Adjust your values

It's hard to say what are sensible values (1.5, 1.1, 1.0001 ?) for any of the ramp-up parameters. It probably strongly differs from case to case and it does depend on the timestep. The recommended approach here is to use the power of empirical research and trial & error.

2. Use modified boundary conditions in b2.boundary.parameters

This is the most complicated step, but at the same time, it might not be important to make the simulation stable (unverified claim). It is, however, important for the simulation to be physically correct. Most of the information here is based on the documentation of the boundary conditions and on the official example ITER_2588_Donly_standalone_drifts from the solps-iter/examples directory. Refer to that example for more details. See e.g. here for a documentation on the boundary conditions. The main idea is that

• There are special versions of the sheath boundary conditions that are modified to properly account for drifts.
• It is advisable to use leakage conditions instead of decay lengths for the radial boundaries. But not sure if that is important for drifts or just a good idea in general.

There are tables below that illustrate transition from "standard" BCs (i.e. those that I happened to be using) to the drift-enabled ones. If you happen to use different BCs, it might be best to check the official ITER example to get an accurate picture of what BCs need to be replaced in your case.

Adjust your values

The values for the BC parameters, such as the leakage factors, are taken directly from the ITER example. They might not be suitable for your case, so make sure to check them and adjust them accordingly.

Boundaries W and E (divertor targets -> sheath): Replace basic sheath conditions by a set of modified ones.

variable	normal	drift-enabled	need to adjust *par
bcene	3	15	no
bceni	3	15	no
bcpot	3	11	no
bccon	3	14	use 1.0
bcmom	3	13	no? (maybe use 0.0?)

Boundaries N (radial from SOL and div region), S (radial from PFRs): Replace decay length conditions by leakage conditions. Leakage conditions are parametrized by a unit-less factor of thermal velocity, i.e. even more arbitrary number. The important part is that it must be negative (and probably smaller than 1). Manual and (ITER) examples seem to be using

• neutrals: -0.01 or -0.001
• ions/electrons: -0.01 for bcene/bcini, -0.001 for bccon

Additionally, in case of potential, make sure to use the zero gradient condition

• i.e. bcpot=2 with potpar=0.0

variable	normal	drift-enabled	need to adjust *par
bcene	9	22	yes (see above)
bceni	9	22	yes (see above)
bccon	9	10	yes (see above)
bcpot	13	2	use 0.0

Boundary S (radial from core)

This one might not be important, but the ITER example uses leakage condition for the core boundary with the continuity equation for neutral particles. The leakage factor used in the ITER example is a bit larger, no idea how to estimate it.

• neutrals: -0.15

variable	normal	drift-enabled	need to adjust *par
bccon	8	10	yes (see above)

3. Use speed-up schema in b2.numerics.parameters

construction This section is under construction.

4. Run and troubleshoot

Under ideal circumstances, the simulation should run smoothly and converge to a solution. More often than not, it will crash, so it is advisable to refer to the troubleshooting checklist below even before you run your first simulation just to save you some trouble. Feel free to revisit it later if (when) you encounter crashes. All in all, there is a lot of moving parts and it is not always clear what could help. Try to prioritize the steps based on the given numbering, quickly go over all of them, and revisit each one later for more thorough treatment. There are some opinionated hints at what seems to be sensible for given step.

Code crash?

Crash in this context means rapid local increase/decrease in some of the plasma parameters and eventual numerical error, which causes the simulation to stop abruptly. In COMPASS simulations, it seems to be often connected to the currents in the divertor region.

1. Start from a converged case without drifts. Speaks for itself, it always helps to start from something else than flat profiles.

2. Start with simpler simulation. Running without Eirene first can help you get acquainted with the process for a smaller computational cost. Same goes for impurities, puffing, and any other complicated stuff. Note that when you enable them later, you might need to redo the whole process from drift factor 0.0. But with a bit of luck, the same settings (timestep, transport, mesh, ...) will work well, letting you to skip the hard part (not verified by experience, yet).

3. Make fluid mesh as even as possible. Check the mesh for any large jumps in cell size and such. If needed, adjust the cell counts and their distribution to fix that and revisit later to do finer adjustments. Consider using finer mesh altogether, make sure to decrease the timestep accordingly in that case. I would not go far over something like 100x50 for a single-null scenario though.

4. Increase the factors gradually. Start with the manual approach, e.g. 0.3, 0.5, 0.7, 1.0. Let the simulation (almost) reach convergence between each step. Try decreasing the increment a bit if it does not work, but no sens in using less then e.g. 0.1 increment. There is also the ramp-up approach, but it seems to be useful only once you know what you are doing.

5. Decrease timestep. Drift simulations typically require smaller dt, but it depends on the speed-up scheme settings (which locally decreases dt by itself). It certainly does not make sense to go much lower than 1e-7, 1e-8. At some point, it actually starts worsening the numerical instability again. And the simulation time becomes impractical of course.

6. Adjust transport. Experience says that drift simulations require modifications to their anomalous transport coefficients. It seems logical, since the drifts are actually one of the contributors to the "anomalous diffusion". But that would seem to dictate a strict decrease in the transport coefficients, which is not always the case. In any case, as you increase the drift factors, keep checking the results and modify the transport coefficients to compensate. Otherwise the plasma parameters might get pushed far from the original state, possibly causing crashes due to "extreme" or "unrealistic" conditions. This is what helped me to reach success in the end, but you can spend an eternity by just looking for the proper non-crashing transport profiles, so it's best to leave any fine tuning as a last resort after you visited the other options. Usually, it's a combination of all of them that helps.