What should we advect with sheared flows

[David Dickinson]

Currently in the ExB shear algorithm in 8.0.6 and 8.1, we have a regular shift of the evolved distribution function g and the potentials phi, apar and bpar in theta0/kx as part of the wave-number remap algorithm.

If we consider that g ~ h - J0 phi then it is possible to note that our current advection approach will change the value of h as both g and phi are kept constant, but theta0/kx changes which will in turn change J0. Similarly, if one takes g and calculates the potentials directly (say using get_init_field) we would expect to find that phi changes after the advection due to the change in kx/theta0.

Should we perhaps simply advect h and then calculate the fields directly from this? The algorithm would be something like

1. Calculate h from g and current potentials
2. Advect h due to ExB
3. Calculate potentials from h
4. Calculate g from h and new potentials.

‌

[David Dickinson]

In ballooning space we shift gnew and the potentials by two pi when the advection takes the distribution function off the end of the grid.

We therefore end up zeroing a two pi section in ballooning space and a full kx cell in a flux tube(1). Currently we zero both gnew and the potentials so this does also set h=0 there.

At the other end of the ballooning space domain we now most likely do not satisfy our parallel boundary conditions for incoming particles. This inconsistency has been observed to introduce grid scale oscillations (in other scenarios, i.e. in runs without flow shear). This may be expected to be damped out after a few time steps, but could potentially cause problems.

(1) If we run a flux tube with nperiod > 1 then this is actually not what we want, instead we just want to zero a 2 pi section of the kx cell.