Exponential decay boundary conditions

[David Dickinson]

Release v8.2.0 (PR #961) introduces an experimental parallel boundary condition, exponential_boundary, which tries to set gnew at the incoming boundary subject to an exponential decay condition. Consider left boundary, the theta point to the right of this is set from gnew(ig+1, 1) = -r * gnew(ig, 1) + ainv * source. If we say

dg/dtheta = - alpha g ~ (g(ig+1)-g(ig))/dtheta then we find g(ig+1) ~ [1 - alpha dtheta] g(ig)

As such gnew(ig+1, 1) = -r * gnew(ig, 1) + ainv * source = [1 - alpha dtheta] gnew(ig, 1)

and hence gnew(ig, 1) = ainv * source / [1 + r - alpha dtheta] provides a way to set the left incoming boundary consistent with an exponential decay.

One must set the exponential decay rate, alpha, and currently this is determined by the user input exponential_boundary_factor. It is likely that the appropriate decay rate will not be a constant as currently implemented and instead will depend on properties of the point in phase-space being considered.

We should consider a smarter way to determine this alpha. One option may be to take the homogeneous GKE (i.e. assume the source terms go to zero). In the collisionless limit this gives dg/dtheta = - [ dg/dt + v_drift . grad g ] / (v|| gradpar). This might suggest alpha ~ [wdrift - i Omega] / (v|| gradpar), Alternatively, writing the time derivative in discrete form we have

alpha g^{j+1} ~ [(g^{j+1} - g^j)/dt + wdrift g^{j+1}]/(v|| gradpar) and hence

alpha ~ [1/dt + wdrift]/(v|| gradpar) - (1/dt) (g^j / g^{j+1})/(v|| gradpar)

Substituting into gnew(ig, 1) = ainv * source / [1 + r - alpha dtheta] gives

gnew(ig, 1)[1 + r - dtheta [1/dt + wdrift]/(v|| gradpar)] + (dtheta/dt) g(ig, 1)/(v|| gradpar) = ainv * source and hence gnew(ig, 1) = [ainv * source - (dtheta/dt) g(ig, 1)/(v|| gradpar)] / [1 + r - dtheta [1/dt + wdrift]/(v|| gradpar)]

Note dtheta [1/dt + wdrift] = (dtheta/dt) [1 + wdrift_func] (roughly).