The linear and collision terms are advanced using the Kotschenreuther algorithm which is essentially a predictor-corrector method. In the time advance steps, operator splitting is used to separate the advances of the linear terms and the collision terms. This is a correct method, but it is not immediately obvious why it works.
The linear terms are advanced using the Kotschenreuther algorithm which is also described in Highcock (2012) and Numata et al. (2010). This is an implicit method, allowing timesteps comparable with ion timescales. The method also avoids inversions of matrices whose dimensions are the problem size, requiring instead inversions of $N_{\vartheta}\times N_{\vartheta}$ sized matrices. This is favourable for parallelization, but is the reason that is kept local to processor.
The algorithm is a predictor-corrector method, with two linear advances to update the
distribution function, and, between these, a field solve to update the electromagnetic
field. These steps are usually called "LFL" (see CMR's note on layouts in
GS2).
Despite this name, the two L
steps are not an operator splitting, as they do not
advance different terms in the gyrokinetic equation. Rather, "(LFL)" is a single operator
which returns the distribution function and fields at the next time step, given their
values at the current timestep.
The electrostatic, collisionless algorithm works as follows (this description is based on Highcock (2012). The discretized electrostatic gyrokinetic equation may be written
where is the distribution function, is the electrostatic potential, and , , and are matrices. The superscript denotes the timestep. For simplicity, we neglect gyroaveraging.
Similarly, we may write the quasineutrality condition as
where again and are matrices. Here represents the velocity space integrals in Maxwell's equations, so contracts with to make a field-sized array.
To obtain the distribution function at the next time step, , we could combine these equations
and solve explicitly for ,
but this involves inverting a matrix of size (where ) at every timestep.
Instead, we use Kotschenreuther's method, which allows us to invert smaller matrices. It is a predictor-corrector method which works as follows. We write
and substitute into the gyrokinetic equation to give,
We then define the predictor and corrector to satisfy parts of this equation independently,
Note that depends only on values at the old timestep, while depends only on the change in the electrostatic potential. The latter can be rewritten as
that is, we can find the response matrix explicitly.
To perform a time step:
We then combine \eqref{eq:defs} and \eqref{eq:response} with the quasineutrality condition \eqref{eq:QN} to obtain the updated field ,
and therefore the field at the new timestep, .
Finally we obtain from the gyrokinetic equation , since we now know , and .
The above three steps are referred to as "LFL". The two linear steps are performed by the function invert_rhs in dist_fn.fpp. These steps are algorithmically the same, since both can be cast as solving the gyrokinetic equation for , given values for , and . That is, the first step is the same as solving the gyrokinetic equation with set equal to . In practice, invert_rhs is always called as:
call invert_rhs (phi, apar, bpar, phinew, aparnew, bparnew, istep)
with g
as a global variable. Note that gnew
is calculated in this function call, but
its input is not used.
Before the first linear step, the fields and distribution function are assigned to the values at the end of the previous timestep:
g = gnew
phi = phinew
apar = aparnew
bpar = bparnew
so that g
contains , and phi
and phinew
both contain , as required
for solving \eqref{eq:gp}.
Before the second linear step, the fields are updated
phinew = phinew + phi
aparnew = aparnew + apar
bparnew = bparnew + bpar
where the field solve getfield
had put the field updates in phinew
,
aparnew
and bparnew
. Therefore phi
now contains , and phinew
contains
. The distribution function g
is not changed before the second linear
step, so g
still contains . Thus calling invert_rhs now solves the
gyrokinetic equation.
These assignments are made in the function advance_implicit in fields_implicit.f90, which then calls timeadv in dist_fn.fpp, which in turn calls invert_rhs.
Note that the corrector is never explicitly evaluated. Note also that the predictor is only used to evaluate the field update in , but is not used again.
g_lo
and gf_lo
memory redistributesThe present default in GS2 is that the linear advances and the field solve are
calculated in the g_lo
layout. This requires no memory redistributes, however the
velocity space integrals in the field solve are very expensive and appear to be the
bottleneck in typical production runs. This led to Adrian Jackson's eCSE project
\cite{JacksoneCSE}, which introduced a new gf_lo
layout for the field solve ,
which improves the field solve at high core counts. This option is also available in the
GS2 trunk. The new approach comes at the cost of a memory redistribute between the first
and , to transform gnew
from g_lo
to gf_lo
.
Note that this only introduces one redistribute, between the first and the
. There is no redistribute between and the second , as the second
uses g
(which is in g_lo
), not gnew
(which is now in gf_lo
).
The above timestepping method is attractive because it only involves inverting and , matrices of size rather than . Introducing collisions would introduce a velocity space dependence into matrix , making it more expensive to invert.
Therefore collisions are implemented using operator splitting. Most of the papers which describe the GS2 algorithm say that collisions are advanced by implicit Euler after doing the whole collisionless time advance (Kotschenreuther et al. 1995 Barnes et al. 2009, Highcock, 2012),
However in the code, collisions are performed between updates in the implicit algorithm, as per Numata (2010) and Roach (2016). That is
To see this in more detail, consider adding a collision term to the discretized gyrokinetic equation, so that
and as before the quasineutrality condition is,
Now inserting the predictor and corrector \eqref{eq:defs} and decomposing the equation as before, we have
It is important to notice that, because is a discretization of a time derivative, it is legitimate to apply operator splitting to equations \eqref{eq:CollisionalGKE}, \eqref{eq:Collisionalgp} and \eqref{eq:Collisionalgc} to advance the collision term separately. Therefore may be obtained as in equations \eqref{eq:L1} and \eqref{eq:C1}, and may be obtained as in \eqref{eq:L2} and \eqref{eq:C2}.
We must also ensure that collisions are included in the field solve. Notice that the algebra in \eqref{eq:FieldUpdate} assumes an expression for the corrector , but once we apply collisions, this expression changes from \eqref{eq:gc} to \eqref{eq:Collisionalgc}. However, this change is automatically accommodated in the code. This is because the rows of the response matrix, , are calculated from \eqref{eq:response} by setting to be a delta function, and noticing that the result of a timestep, , is the required row. Therefore, all the terms contained in the timestep are automatically included in the response matrix.
For the field solve to be consistent, it is therefore only necessary for the timestep used
in calculating the response matrix to be the same as the timestep used in calculating the
predictor, . That is, in timeadv, we must call the same routines for
istep == 0
(the response matrix) as for istep /= istep_last
(the predictor step).
Introducing the collision term introduces four new memory redistributes, in addition to the one redistribute in the collisionless case. In all, we have five redistributes at the following points:
gnew
from g_lo
to le_lo
to calculate collisions on
the predictor .gnew
back from le_lo
to g_lo
after collisions.gnew
from g_lo
to gf_lo
to use the predictor in the
field solve; there is no direct redistribute from le_lo
to
gf_lo
.gnew
from g_lo
to le_lo
to calculate collisions on
.gnew
back from le_lo
to g_lo
after collisions.In the above discussion we have neglected the fact that collisions are applied on \(h\) rather than \(g\). This does not affect the operator splitting in the Kotschenreuther algorithm, but it does affect the splitting between the collisionless and collisional advance.
A problem arises with g_adjust if collisions are fully separated from the collisionless Kotschenreuther advance. The method may be summarized as
That is, collisionless versions of the new distribution function and fields are found from a collisionless Kotschenreuther advance, then g_adjust is used to find a collisionless , collisions are applied, g_adjust applied again to find the new distribution function , and then finally a field solve is called to find new fields consistent with the new distribution function.
This is the current approach in the "split collisions" version of GS2. However, it is not quite correct, as the second g_adjust uses fields which are consistent with not . The simplest correction to this would be to find the new fields with a field solve on before calling the second g_adjust,
This method requires a new field solve to be written for , but like which is used in the initialization, it would not involve the response matrix.
The new field solve would also only need to be called in electromagnetic runs, as and are (probably) the same by conservation of mass.