Error control for the split nonlinear scheme

[David Dickinson]

With the split nonlinear approach due for inclusion in 8.2 we can separate the evolution of g due to the nonlinear and linear operators. This operator splitting can be helpful but introduces an error ~ O(dt). Currently we don’t have anyway to track this error. Ideally one could calculate this and adapt the linear time step to keep the error below some threshold. This issue is to discuss approaches to this.

‌

Firstly let us schematically write out the original algorithm:

dg /dt = L + NL --> g_{i+1} = g_{i} + dt * L(g_{i}, chi_{i}, chi_{i+1}) + dt * N(g_{i}, chi_{i})
with L the linear operator, N the nonlinear operator, chi the gyrokinetic potential and subscripts represent the time index

In the split scheme we first evolve g using only the nonlinear operator:

g*_{i+1} = g_{i} + dt * N(g_{i}, chi_{i}) (Note strictly the nonlinear operator may be advanced with many sub-steps)

We then evolve the linear term

g'_{i+1} = g*_{i+1} + dt * L(g*_{i+1}, chi_{i}, chi'_{i+1})

where we note that strictly chi_{i+1} /= chi'_{i+1}

‌

If we consider the unsplit scheme to be the true answer, we can estimate the error, e, as

e ~ g_{i+1} - g'_{i+1} = g_{i} + dt * L(g_{i}, chi_{i}, chi_{i+1}) + dt * N(g_{i}, chi_{i}) - g_{i} - dt * N(g_{i}, chi_{i}) - dt * L(g*_{i+1}, chi_{i}, chi'_{i+1}) = dt * L(g_{i}, chi_{i}, chi_{i+1}) - dt * L(g*_{i+1}, chi_{i}, chi'_{i+1})

where we have assumed that the update from the nonlinear term is the same in both split and unsplit cases (not strictly true).

We can write the intermediate g*_{i+1} as g*_{i+1} = g_{i} + dg and using the linearity of L we get

e ~ dt * L(g_{i}, chi_{i}, chi_{i+1}) - dt * L(g_{i}, chi_{i}, chi'_{i+1}) - dt * L(dg, chi_{i}, chi'_{i+1})

If one assumes chi_{i+1} = chi'_{i+1} (again not true but perhaps not unreasonable when NL term significant) then the first two terms cancel and we’re left with

e ~ - dt * L(dg, chi_{i}, chi'_{i+1})

that is, the error is approximately the change in g we get from taking one full linear time step starting with g = dg. If one assumes that the change in chi from the field solve at the mid step of the full linear time step is small (i.e chi'_{i+1} ~ chi_{i}) then one finds e ~ - dt * L(dg, chi_{i}, chi_{i}). It is then possible to calculate this error estimate by setting g = dg_nonlinear , calling invert_rhs to calculate the corresponding update to g such that error ~ gnew - g.

This error estimate requires one additional call to invert_rhs and an additional distribution function sized array to store dg (transiently). As such it should be relatively cheap. This approach is built on many assumptions/approximations and also entirely neglects any potential effect from collisions (to improve this one may replace the call to invert_rhs with effectively a call to timeadv without recalculating the NL terms).