Dubious response from A||

[David Dickinson]

Setting up a ballooning space simulation with only A|| one can find a strange behaviour in the response matrices. Here is the raw response (i.e. pre-inversion) for a case with nperiod = 1, ntheta = 32 and delt = 1e-5.

Moving in the horizontal direction corresponds to considering different locations where A|| is perturbed, whilst moving in the vertical direction corresponds to considering the response at different theta locations for a fixed location where A|| is perturbed. One can see two vertical stripes, one at each boundary. These correspond to the response at all theta when A|| is perturbed at the boundaries. These boundary stripes remain for arbitrarily large nperiod and also in simulations with other fields included.

The response here is fieldeqa = antota - kperp2 * aparnew and antota ~ Sum_species(q_s * Int (sigma * |v_par| * gnew * J0)). The data in the plot is mostly dominated by the antota contribution. The diagonal outside of the boundaries is ~370, whilst kperp2*A|| ~ 4e-3, off diagonal antota is close to 1e-2. At the boundaries the response in the vertical stripes is all around 185.

Consider a perturbation in aparnew at ig = -ntgrid which means the source is only non-zero at ig = -ntgrid. We first integrate for gnew with sigma = -1 from ig = ntgrid to ig = -ntgrid, starting with gnew = 0. Clearly this gives `gnew = 0` for the range ig = [-ntgrid+1:ntgrid] as the source is zero. At the last point we can get a non-zero gnew value. Next we integrate from ig=-ntgrid up to ig=ntgrid for sigma = +1. Here we again start with gnew = 0, but encounter a non-zero source immediately and therefore can find gnew /= 0 for the full theta range (except right at the boundary). This difference means that it’s not possible for the two signs of v_par to cancel for most theta here, which could explain these vertical stripes.

When we perturb away from the boundary, say ig = -ntgrid + 1 , we might expect a similar behaviour but we’re not obviously seeing this here. All of the values above the diagonal must be due to gnew for sigma = +1 but there is a significant difference for perturbing on the boundary or next to it. This may suggest that there’s some substantive difference in the algorithm between these two situations, although it’s not obvious what this may be.

‌

‌

[David Dickinson]

@Colin Malcolm Roach may be interested in this.

I originally thought that this might be something to do with non-zero currents at the parallel boundaries, but this behaviour remains even if I construct a zero boundary current solution using the homogeneous solutions.

[David Dickinson]

Here the relevant source is dA||/dt. We calculate this as:

‌

        apar_m = (aparnew(ig+1,it,ik) - apar(ig+1,it,ik)) + &
             (aparnew(ig,it,ik) - apar(ig,it,ik))

‌

This shows how we calculate the source at the grid centres. Note how for a boundary perturbation this source is only non-zero at a single theta grid point (the boundary), whilst for all other places we get the source appearing at two theta grid locations. If we arbitrarily remove this and only evaluate the source on the grid the response is very different:

‌

‌

Excluding the last column and inverting we end up with something fairly strongly diagonal. The last column is now dubious as this corresponds to the response from perturbing the last theta grid point, but the source term is not currently evaluated here so it only gets into the system usually due to the grid centring. Ignoring the last row and column we get a response which looks like

‌

Perhaps we should consider introducing a ghost point for our boundaries rather than enforcing them on the grid which is evolved? This would be equivalent to having a fictitious point at ig = -ntgrid - 1 where we enforce our boundary conditions (gnew = phinew = aparnew = bparnew = 0.0) and where we start our parallel integration. We would not consider perturbing the fields at these fictitious points of course.

[David Dickinson]

Increasing the upwinding (bakdif = 0.5) seems to give a better response

‌

But this actually has the same issue when inverted:

(whilst inverting excluding the first and last rows/columns gives something strongly diagonal).

‌

[David Dickinson]

Some motivation for justifying why we might care about this (beyond it looking dubious) is from the condition number. For this test case the condition number of the response matrix is about 1.5e6 whilst the condition number for the sub-matrix dropping the first and last rows is 1.004

[David Dickinson]

The vertical stripes are closely associated with the r*gnew part of the rhs.

The equation we’re integrating in invert_rhs_1 is effectively dg/dtheta = Source which when discretised gives something like g(ig+1) = source(ig) - g(ig) . At the boundary g(ig) is forced to zero so we have g(ig+1) = source(ig) at the next point we get g(ig+2) = source(ig+1) - g(ig+1). When the source is non-zero away from the boundary (i.e. we’re perturbing A|| away from the boundary) then the grid centring gives us the same source at adjacent points. This means that we’d have something like g(1) = source(0) - g(0) = source(0) and then g(2) = source(1) - g(1) = source(0) - g(1) = source(0) - source(0) = 0. In other words, this works out such that we only get a response where the field is perturbed and it is zero everywhere else. This would avoid vertical bands. (Of course this is an idealised situation and we would not end up cancelling exactly, rather just mostly cancelling), At the boundary we end up with a different situation; g(ib+1) = source(ib) - g(ib) = source(ib) and g(ib+2) = source(ib+2) - g(ib+1) = 0 -g(ib+1) = 0-source(ib) and we’d end up with something effectively oscillating in theta, which is what we see (the previous plots all take the magnitude of the response so it can be put on a log scale, if we don’t do this then it alternates sign on the grid scale).

‌

If we had a fictitious point at which we set our boundary condition then g(ib) would be non-zero and this should look like a non-boundary point (i.e. we end up mostly cancelling the response away from this point). A quick attempt to hack this in does avoid the vertical stripes at the boundaries, but introduces them at the points just in from the boundary. This is not too surprising as I didn’t extend the source, I just mirrored the boundary sources.

[David Dickinson]

Hacking in something like a fictitious boundary point in order to set the gnew values at the incoming theta grid point gives a response matrix which looks like

‌

and has a condition number very close to one.

This effectively assumes that we have gnew = phi = apar = bpar = source = 0 at -ntgrid-1 and ntgrid+1, calculates the effective source between this fictitious point and the first/last theta grid point and then uses this to set gnew at the start/end of the grid. This may be harder to extend to the nonad_zero case where set hnew = 0 but allow phi and bpar to be non-zero. This would require us to somehow find phi and bpar on this fictitious point. This could perhaps be done via extrapolation.

‌

Ultimately I think this issue stems from the fact that we’re currently perturbing the potentials at the same location as we enforce the boundary condition and thereby constraining how the distribution function can respond. This then can lead to significant differences in the response structure at the boundary.

‌

[David Dickinson]

Whilst the current implementation is quite hacky (and probably not correct) it looks like it could possibly improve the behaviour of the code in some cases, for example here are comparisons of omega and gamma vs time for a simulation with only A|| . Note how much smoother the time traces are are the hacked version.

‌