Drifts

Feature blog entries introduce advanced concepts

The reader is assumed to know how to create a simple SOLPS-ITER simulation and to possess some user wisdom.

Quick checklist to enable drift terms:

Enable drifts in b2mn.dat using the b2news_fac* switches.
Modify boundary conditions in b2.boundary.parameters to account for drifts.
Apply the speed-up scheme in b2.numerics.parameters to improve numerical stability and computational speed.
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

Charged particles in non-uniform magnetic and electric fields are subject to drifts of the guiding centre. For SOLPS-ITER modelling, the most important drifts are the diamagnetic drift and the ExB drift.

On the diamagnetic drift nomenclature

In SOLPS-ITER literature and documentation, you'll encounter the diamagnetic drift under three names: the diamagnetic drift, the grad-B drift and the Pfirsch-Schlütter currents. Here's why.

The Braginskii equations, which lie at the heart of the B2.5 plasma solver, have the form time derivative + divergence of flux = sources. Drift terms appear mostly inside the divergences, as means of transporting particles, momentum and energy from cell to cell. Now, it can be mathematically shown that the diamagnetic drift can be separated into a part which has zero divergence (is "divergence-free") and non-zero divergence (which is essentially the grad-B drift). While implementing drifts in B2, it was found that using the entire diamagnetic drift inside the divergence leads to numerical instabilities. So a trick was performed: the divergence-free part was excluded from the equations. This made the drift simulations less likely to crash and burn, but it isn't the whole story yet.

In a tokamak, the grad-B drift is vertical and of opposite direction for electrons and ions, so it causes vertical charge separation. The resulting electric field drives electric currents along the magnetic field lines. These are the Pfirsch-Schlütter currents, and their divergence is the same (taken negatively) as the divergence of the grad-B drift. These currents are even easier to calculate than the grad-B drift inside the B2.5 equations. So a final transformation was made where the originally diamagnetic term was replaced by "Pfirsh-Schlütter flows" inside the divergences.

Agonisingly, the physical nature of these three phenomena is all different. Diamagnetic flows are cross-field and fluid, the grad-B drift is cross-field and for the guiding centre, and P-Sch currents are parallel. If you calculate them on their own, their magnitude and direction is all different. The only thing they have in common is that their divergence is the same. So, keep in mind: inside divergences, diamagnetic, grad-B and Pfirsch-Schlütter are interchangeable, but if you're calculating real particle fluxes (like for energy convection), you have to use the entire diamagnetic drift formula.

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:

Constant profile \(f(t)=const.\)
set _inc = _dec = 1.0, then \(f_0 = f_n = f_{target}\) is kept during the whole run
increase _start = _target = 0.3 manually between runs (e.g. 0.3, 0.5, 0.7, 1.0)
Ramp-up profile, customize how steep
set _start = 0.1 (something small), _target = 1.0
set _inc = 1.5 (something larger than 1)
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.

Start from a converged case without drifts. Speaks for itself, it always helps to start from something else than flat profiles.
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).
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.
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.
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.
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.