Linearly unstable zonal modes

[David Dickinson]

The zonal modes should be linearly stable, however in a number of simple test cases this isn’t true and the zonal modes are subject to a numerical instability.

‌

Attached is a simple input file for testing/demonstrating this.

‌

The base case is a two field (phi+apar) sheared cylinder with no grid decentring. The figure below demonstrates the numerical instability and shows that the growth rate is sensitive to the inverse aspect ratio and in the slab limit (when there are no trapped particles) this instability goes away.

‌

The instability is essentially a parallel grid-scale mode. This should be suppressed by some grid decentring and the below figure illustrates that this is indeed effective in some cases. In particular with opt_source = .false. and bakdif = 0.05 for the two field case of for the phi only case with bakdif = 0.05 and either opt_source the mode seems to be suppressed. Note that opt_source = .false. means we have spatial de-centring on the apar source terms as well, however this has been seen to cause problems linearly for non-zonal modes as discussed in issue #91 and issue #97.

Another option is to keep bakdif = 0.0 but introduce time decentring and this is shown below. It can be seen that this does help stabilise the numerical instability and the required value of fexpr is sensitive to the timestep. This highlights a potential issue for nonlinear simulations. One may do linear analysis and pick a fexpr that’s sufficient to prevent this numerical mode, however as the amplitude grows and the time step falls in the nonlinear simulation, the effective stabilisation drops and this numerical instability could start to grow. This will cause further time step reductions and related reductions in the stabilisation resulting in a feedback cycle leading to the simulation blowing up with a tiny time step.

This probably impacts a range of simulation types, but I think is likely to be more of an issue in EM cases as we can’t rely on bakdif to provide some time independent dissipation without having an impact on the linear physics.

‌

For this simple problem the field source terms are essentially just v_par grad_par(phi) + v_par dA/dt, i.e. v_par E_par and the system is just a periodic slab/cylinder so the problem is fairly simple.

‌

[David Dickinson]

Looking at the response matrix for the two field, bakdif = 0.05, eps=0.01 case pre-inversion there are certain grid points where the response is very close to zero. This is likely due to symmetry in the response etc. Setting eps = 0 enhances this further. These points arise in the apar-phi and apar-bpar mixed parts of the response (i.e. off-diagonal blocks).

There are also zig-zag like properties in the response matrix in the diagonal blocks (e.g. in the phi-phi block). This gets better with increasing ntheta, but remains near the boundary even for quite high ntheta. This could perhaps exacerbate/seed zig-zag instabilities appearing in the fields.

‌

[David Dickinson]

This could be related to issue #146 where the velocity space errors can become large when bmag is close to bmax.

[Michael Hardman]

Is this related to the numerical instability featured in appendix A of https://doi.org/10.1016/j.jcp.2019.01.025 ?

Michael Barnes was able to eliminate a grid scale numerical instability by using a different differencing for the distribution function as for the source.

[Michael Hardman]

I imagine that this could be tested for "range" geometry relatively simply, although generalisation to "box" geometry might be more challenging.

[David Dickinson]

It is possible that this is related to the case described in the stella paper, although here we only seem to see the instability when including apar alongside phi (i.e. not in electrostatic simulations). I think this is also likely to be related to issue #142 which demonstrates how an instability can only appear once we go below a critical time step.

[Michael Hardman]

What does the green curve in the second figure above represent? Before your comment, I thought it indicated that the numerical instability grew with only phi =/= 0.

To see whether this apar issue was related to the large size of v|| for electrons, I ran with an electron mass =1. The numerical instability remains. If d apar / dt is the problem, would solving for g = h - <X> Z e F_0 /T = < delta f> rather than g_gs2 be interesting?

‌

[David Dickinson]

Ah yes, sorry I seem to have misremembered this – you’re correct that the green curve in the second figure suggests that the instability can persist if phi is the only field included. I will have to check if I can reproduce that case, it may be an incorrect label.

‌

I think that would certainly be interesting. I have some concerns about our treatment of d apar / dt currently (this is forced to zero for the first linear solve of each time step, for example).

[Michael Hardman]

After checking my output for longer times, it seems fapar = 0reduces the growth rate of the numerical instability by several orders of magnitude, but does not kill it entirely. A numerical instability appears to exist just with phi. If I have time I will try the differencing that Barnes used in stella, and look at the apar term.

[David Dickinson]

Zeroing the trapped pitch angle integration weights seems to avoid the blow up, helping further suggest that this is at least partially related to the trapped particles. Zeroing the trapped particle source in invert_rhs_1 delays but does not avoid the blow up and grid scale instability entirely. Removing the source will leave gnew as 0 for trapped particles after calling invert_rhs_1, however when we operate on h (i.e. when we call g_adjust) this will be non-zero so this isn’t the same as removing the trapped contribution to integrals if we were to operate on h . For the collisionless linear problem I don’t think we would expect to be calling g_adjust for anything, but may be worth exploring.

Increasing ntheta significantly seems to help delay the issue slightly.

[Michael Hardman]

That is interesting. My experience with forcing the passing electron response to be zero for the ‘modified adiabiatic response’ option would suggest that setting h = 0 is not the same as removing h from the problem. In order to get a robust solution I had to also enforce h = 0 properly in the response matrix calculation.

[David Dickinson]

Yes I might try this approach as well.

‌

I also note that we can delay the onset of the issue by increasing the velocity space resolution somewhat. Below are two cases where I increase npassing and nesub significantly from the baseline case also shown (blue).

‌

[Michael Hardman]

I am surprised that changing npassing has an effect, given that the problem appears to be associated with trapped particles. When you changed ntheta, did you see a similar delay in the onset of the instability. It is possible to separate the effects of increasing ntheta from increasing the # of trapped pitch angles? I recall that Jason Parisi implemented such a feature, although it was abandoned because nonlinear simulations tended to go numerically unstable if too few trapped particle pitch angles were included.

[Michael Hardman]

How are the trapped particle integration weights tested currently? Is it possible that the problem arises because of the trapped particle contributions to the velocity integrals of the Maxwellian in the field solve? This seems unlikely, but perhaps consistent with what you have observed.

[David Dickinson]

Yes I this may be due to the error on the field integration which then feeds back on time advance for all pitch angles. Perhaps the trapped algorithm is particularly sensitive to this for some reason.

We currently test the velocity integration by integrating a few functions with known results, although this tests the weights across the entire velocity space. We also check that the sum of the weights in different regions (integrating one) and report the error in <run_name>.vspace_integration_error

‌

Increasing ntheta did similarly delay the onset.

[David Dickinson]

Comparing the fields at the end of the time step to those calculated directly from gnew (get_init_field) shows a relative error of ~1e-7 suggesting that there’s no major inconsistency in the response matrices (unsurprising result).

[David Dickinson]

I wonder if this is related to issue #144. There (and possibly in other issues) it can be shown that, at least in the case with no drifts (as here), the homogeneous solution under our current assumptions and with our space and time grid-centre approach can lead to grid scale instability. In particular one can find a CFL like condition where dt v|| / dtheta > 1 to avoid a change in sign between theta grid points. This is actually the opposite of a typical CFL limit in that it requires large distances from the advection compared to grid spacing, however this is because we are effectively integrating in theta and not time (so dtheta plays the role of dt in a typical cfl condition).

‌

This condition seems to match reasonably well with some of our observations:

1. The instability is delayed/suppressed with larger dt.
2. The instability is delayed/suppressed with larger ntheta (smaller dtheta).
3. The instability appears to require trapped particles – these particles experience low v|| (v||=0 at the bounce point) so are likely more sensitive to this.

My working hypothesis is that the homogeneous solution for trapped particles contains grid scale at least close to the bounce point. This then provides a seed for grid scale in the potentials which can then infect all of velocity space, reinforcing the issue and leading to a numerical instability with high k||. I’m not sure that this can explain all of our observations yet (e.g. dependence on velocity space resolution etc.).

One assumption that we appear to make in the homogenous solve is that g_h at the previous time step is exactly zero. This could potentially cause problems as it doesn’t actually satisfy the boundary conditions that we impose (g_h=1 at the bounce point) but setting g_h at the previous time step to exactly one doesn’t change the current behaviour.

The only way within the current approach to guarantee that the grid scale will be avoided for all v|| / all pitch angles is to set bakdif = 1. The solution for g_h then is 1 at the bounce point and zero everywhere else (i.e. over damped). This is clearly not a useful solution as we cannot then use g_h to ensure continuity at the bounce points in the final solution.

It appears we may need to treat the homogeneous solution a little differently than the regular solve, particularly near bounce points. It may be sufficient to adopt special handling near bounce points, but one can imagine that for sufficiently low energy particles one may end up violating the CFL-like condition across a broad range of theta.

[David Dickinson]

Some counter-evidence – if I scale the parallel velocity down (using vparknob) then the instability occurs later rather than earlier despite the condition in the previous comment suggesting this makes the grid scale more likely. Of course v|| enters in more places than just this (and for example scales the A|| source term which we know is also important for this blow up).

[David Dickinson]

To help with this discussion I’ve written out the discretisation.

‌

We are considering dg/dt + v|| grad|| g = 0

We write dg/dt = (1/(2dt)) [(1+b) [g_{i+1}^{n+1} - g_{i+1}^n ] + (1-b) [g_{i}^{n+1} - g_{i}^n ] where g is our homogeneous solution, b is bakdif, subscripts represent theta grid points and superscripts represent time indices. This is dg/dt evaluated (approximated) at the theta grid ~centre, upwinded by b.

We write grad|| g = dg/dtheta = (1/dtheta) * [ (1-f) [g_{i+1}^{n+1} - g_{i}^{n+1}] + f [g_{i+1}^{n}-g_{i}^{n}]] (note for ease we take gradpar=1, in general there will be a gradpar = dtheta/dl factor included in the coefficient).

‌

Using these expressions in our equation and defining s = 2 dt v|| / dtheta and collecting terms we get

g_{i+1}^{n+1} [1+b + s(1-f)] +
g_{i}^{n+1} [1-b - s(1-f)] +
g_{i+1}^{n} [-(1+b) + s f ] +
g_{i}^{n} [-(1-b) - s f] = 0

‌

If we let A = [1+b + s(1-f)] then we can write the above as

g_{i+1}^{n+1} A +
g_{i}^{n+1} (2-A) +
g_{i+1}^{n} (s-A) +
g_{i}^{n} (A - s  - 2) = 0

Dividing through by A and writing as an equation for g_{i+1}^{n+1} we get

g_{i+1}^{n+1} =
g_{i}^{n+1} (1-2/A) +
g_{i+1}^{n} (1-s/A) +
g_{i}^{n} ([s+2]/A - 1)

‌

Now if we assume that g^n = 0 (as done in the code, also the same result if we simply assume g^n==const we get

g_{i+1}^{n+1} =
g_{i}^{n+1} (1-2/A) +

One can see that if 2/A > 1 then we end up with g at the next theta grid point being the opposite sign to the current grid point, a clear sign of grid scale (if 2/A >1 remains true for all theta then we get grid scale which may be damped or growing depending on if 2/A > 2 or not).

Taking bakdif = 0 and fexpr = 0.5 (our second order scheme) we end up with A = 1 + s/2. When s->0 we get A=1 and hence g_{i+1}^{n+1} = -g_{i}^{n+1} – an exact grid scale response. If we make bakdif=1 then A=2 when s->0 and we end up with g_{i+1} ^{n+1} = 0 – over damped to zero (and not useful for the intended purpose).

Of course in our usual treatment we do not use v|| on the grid but the grid centre values so we shouldn’t get exactly v|| (and hence s==0) in the numerical implementation. If we write s_min = (2dt/dtheta) v||_min where v||_min is the minimum grid centre v|| for the pitch angle of interest then we can find the following requirement to avoid grid scale:

2/A < 1 --> 2/[1+b+s_min (1-f)] < 1 -> 2<[1+b+s_min(1-f)] = 1 + b + (1-f) (2 dt / dtheta) v||_min] -> (1 - b) * (dtheta / [2 dt (1-f)]) < v||_min

Taking fexpr = 0.5 this simplifies to (1-b) * (dtheta / dt) < v||_min. As bakdif if a constant for all pitch angles we need to choose b such that 1 - v||_min dt/dtheta<b<1

In the code this would look something like 1 - min(vpar / spec%zstm) < bakdif < 1 (where the minimum ignores forbidden regions)

‌

Next steps

The first obvious step is to check the above condition and report when it isn’t satisfied to see if this is linked to when we see this instability. We may be able to tolerate a little bit of grid scale in the homogeneous solution as we may only add a small amount of this into the final solution, trapped particles may not be significant etc. As such this condition being violated may not mean we see grid scale arising.

Following that one may want to try to find a way to avoid this limitation or at least reduce the severity. One avenue is to remove the assumption that g_n is constant, although it’s not clear that this will help. Another avenue is to explore the use of a different approach to centering.

At least in this simplified case one may be able to assume an eigenmode (d/dt → -i omega) and use an analytic solution (e.g. https://www.wolframalpha.com/input?i=Solve[-i*omega*f(x)+%2B+v(x)+d+f(x)%2Fdx+%3D+0]) although care would still be needed near the bounce point. Again, it’s not clear that this is viable/appropriate in general.

Ultimately we are using an explicit scheme to integrate in theta and this will be subject to CFL like restrictions. Can we switch to an impliciit scheme?

Godunov’s theorem also tells us that only first order schemes avoid the introduction of new turning points in advection problems (assuming a linear scheme). Our grid centre approach to the parallel derivative here means we end up with essentially the second order accurate central difference parallel derivative evaluated at the grid centre (when b = 0).

[David Dickinson]

For this test case I find we would need bakdif > 0.99997035125329015 to satisfy the above condition (min v|| ~ 1.0e-4). This does indeed avoid the numerical instability but one may assume that this is going to be excessively dissipative.

Perhaps using large (~1) bakdif just for the homogeneous problem may allow us to avoid the grid scale in the homogeneous solve whilst not being overly dissipative in the full solve?

[Michael Hardman]

Thank you for writing out these notes carefully. It seems that you might have discovered the underlying issue at last! I had not noticed the assumption that g^n = 0, so I shall have to reread these code blocks.

If your suggestion to force bakdiff = 1 in specific places in the code works, then I think that we would need to do a convergence study of some kind to see if the code is noticeably first-order accurate as a result. If you were able to show that eigenfunctions and growth rates converged rapidly enough, then I think I would be happy to go forward with this, especially if it allowed GS2 to be used for some of the cases where we have previously had issues.

[David Dickinson]

Yes we would definitely need to do some careful testing of this, but I am hopeful that it shouldn’t impact the convergence too much.

I’ll have a go at implementing this today and will report back here.

[David Dickinson]

Good and bad news. On the bad side adding this change seems to increase the growth rate of the instability:

On the positive side it does seem that with this change we are avoiding pure grid scale taking over:

This was quite a quick hack implementation so it’s possible I’ve missed something that is now inconsistent (although I think it should broadly be correct).

[Michael Hardman]

Having reflected on what I remember of the parallel solve, I wonder if it is helpful to consider the origin of the `homogeneous' problem. I think that I have convinced myself that setting g^n = 0 is in fact the right thing to do in the parallel solve of the `homogeneous' part of g.

Overall in the parallel solve, my understanding is that we need to integrate an equation of the form

d g / dt + v|| d g / d theta = S

where S = S(t,theta,phi,g) is some source term. We need to ensure continuity of the solution along theta across the boundaries between data that is stored in different parts in memory, as well as ensure that periodic BCs are satisfied for passing particles on zonal modes, and bounce conditions are satisfied for the trapped particles. To overcome this problem we split the problem into a `homogeneous' and `inhomogeneous' part so that we can treat the sources and the boundary conditions as separate problems.

I find it helpful to first discretise in time, and then argue for how to split the discretised equation.

g^n+1 - g^n + Dt v|| [(1 - f ) d g^n+1 / d theta + f d g^n / d theta ] = (1 - f ) S^n+1 + f S^n

Assuming that phi^n+1 is known via response matrix methods, (and ignoring the drifts in S^n+1) the above equation is a straightforward ODE in theta for g^n+1, that we can rewrite as

g^n+1 + Dt v|| (1 - f ) d g^n+1 / d theta = RHS(known fields)

We choose to split g^n+1 into a inhomogeneous bit g^n+1_IH and a homogeneous bit g^n+1_H. One part has the known source, the other has a nonzero incoming boundary condition.

g^n+1_H satisfies

g^n+1_H + Dt v|| (1 - f ) d g^n+1_H / d theta = 0 (with g^n+1_H = 1 at theta = Lmin)

g^n+1_IH satisfies

g^n+1_IH + Dt v|| (1 - f ) d g^n+1_IH / d theta = RHS(known fields) (with g^n+1_IH = 0 at theta = Lmin)

Once we have the two solutions g^n+1_IH and g^n+1_H, we add them together in the right ratio to satisfy our boundary conditions.

Additional comments: @David Dickinson what do you mean by ‘explicit scheme to integrate in theta’. Isn’t the scheme implicit by virtue of having points evaluated at theta_i+1 as well as theta_i in both terms in the LHS of the equation? Perhaps I am missing the point.

I like your explanation of the grid scale oscillation above, as I think it explains why it is the trapped particles that have appeared to be problematic, even though all particles nominally have the same parallel integrator. Trapped particles are the only ones whose velocity → 0 on the grid! However, this problem would appear to get worse with increasing ntheta resolution, as that would put the minimum v|| used by the integrator closer and closer to zero at the bounce point. Is this consistent with your observations?

[David Dickinson]

I think this works out to be an explicit scheme as, once we have dropped g^n, we are solving for g^{n+1}_g^{i+1} using only fixed quantities on the grid and g^{n+1}_{i}.

‌

The problem appears to get better for this test case with increasing resolution. Note that the term which goes to 0 as v|| → 0 contains v||/dtheta, so it may depend which goes to zero faster as to what behaviour we expect to find.

[Michael Hardman]

I would agree that the scheme is explicit in theta if we took bakdiff = -1 in the parallel integrator (keeping only i+1 points from the d/dtheta). Do you mean explicit in time?

Does your final observation give a route to deal with the instability? If we choose sqrt{1 - lambda B(theta_bounce +- D theta)}/ D theta appropriately, perhaps we can guarantee that v||_min is bounded from below. If this was the problem, this would only work for a fixed energy grid since v|| ~ sqrt{energy}.

[David Dickinson]

Ah yes I see your point about explicit vs implicit here, sorry!

We might be able to do something with careful choice of grid points, however this is likely to be complicated as the bounce point for one pitch angle is the point one grid point away from the bounce point of the previous pitch angle, so the placement would have to be appropriate for all pitch angles simulataneously.

‌

One issue with the modified approach of full upwinding for g_h is that if this gives strong damping then g_h will be close to zero at the bounce point on the other side of the orbit from where we set g_h=1 and start integrating. This can mean that we need to add in a rather large amount of the homogeneous solution into the full solution, resulting in a large dg_trapped/dt due to this.

[David Dickinson]

I’ve added a routine to calculate the homogeneous solution and write this in a restart file format. This can be used to verify that for bakdif << 1 the homogeneous solution changes sign each grid point. For the test case at hand this remains true until bakdif ~ 0.998. By looking at the magnitude of g_h we can see that we lose much of the structure in the solution when bakdif is large – we get pure exponential damping (straight line on the log plot) with a large damping coefficient.

‌

[David Dickinson]

It seems we can suppress the grid scale in the solution using this treatment, but we can’t avoid numerical instability. The instability still appears due to the trapped particles (I have only implemented the bakdif change for trapped particles and zeroing the trapped particle integration weights again sees the instability go away in phi2).

If I also force bakdif = 1 for the inhomogeneous trapped particles then we do suppress the growth quite a lot but grid scale reappears.

‌

I’ve pushed experimental/upwind_trapped_homogeneous which includes the change to the trapped homogeneous solve and the routine for saving the solution.

‌

[David Dickinson]

I’ve also been looking at the full solution at the end of the simulation. I noticed that the gnew tends to have a lot more theta structure than the “g_wesson” piece, h. There’s also an interesting and obvious asymmetry in the mode structure / dist-fn even thought the problem should be entirely symmetric.

‌

[David Dickinson]

Given we are not really advancing the homogeneous solution in time and we can add in an arbitrary amount of the resulting solution, is there any reason we can’t choose dt arbitrarily large when solving the homogeneous problem? For sufficiently large dt we should be able to ensure that the required condition is satisfied everywhere.

[David Dickinson]

The source term used for the inhomogeneous problem is made from the actual RHS of the GKE plus the terms coming from the advection of the current g. Interestingly, zeroing the terms from the RHS doesn’t avoid the numerical instability whilst zeroing the advection terms (just for trapped particles) does avoid instability provided a small amount of bakdif is used. These advection terms (a*g+b*g) look a lot like terms appearing in the homogeneous problem and one may therefore anticipate that they might provide another possible source of grid scale etc. On their own, these terms shouldn’t lead to growth.

[David Dickinson]

One may note that for this problem the source terms arising from the potentials are all multiplied by v|| so will approach zero near the bounce points. The inhomogeneous problem then resembles the homogeneous one (with the same ability to generate grid scale) with the primary differences being that the inhomogeneous problem includes the g from the current time step (which we currently take to be zero in the homogeneous treatment) and we integrate in theta with different initial conditions (we start with gnew = 0 for the inhomogeneous problem and g1 = 1 for the homogeneous one).

[David Dickinson]

Unfortunately I managed to introduce a bug whilst adding some tests into the code so results shown since the entry regarding the experimental branch should be ignored. The experimental branch does not contain the bug.

‌

To summarise, what I think I’ve learnt is

1. The trapped homogeneous solution is almost guaranteed to be a (lightly damped) grid scale near bounce points
2. Zeroing the integration weights for trapped particles avoids zonal phi growing in time
3. Increasing the time step seems to avoid the issue
4. Instability remains if one swaps ky and kx (i.e. not exclusively an issue with zonal modes -- not due to periodic b.c.)
5. Instability can be somewhat suppressed with bakdif, but even by bakdif = 0.5 the mode is still weakly growing. With bakdif = 0.95 it decays slightly (I've not narrowed down the transition any further). This is with opt_source = T
6. Turning off A|| avoids the growth provided bakdif > 0 (with bakdif = 0.0005 it is stable for a while but then blows up at t~5, bakdif = 0 blows up at t~1, bakdif = 0.005 is fine all the way to t=100 at the end of the simulation -- some signs it is perhaps about to go bad).
7. Keeping A|| and bakdif = 0.005, we can suppress the instability by switching to opt_source = F, which adds bakdif dissipation onto the A|| source term. For this problem it is simply adding bakdif to the dA||/dt term.

Point 1 is possibly a red-herring. This is pretty much always true but doesn't always cause problems. Not impossible that this can cause some problems but we could avoid introducing gridscale if:

1. The inhomogeneous solution doesn't need much homogeneous solution as a correction to fix the boundary conditions (probably true for cases with weak source)
2. Trapped particles don't contribute much to the potentials etc.

Of course, one might expect the inhomogeneous solve with a weak source to resemble the homogeneous problem and hence we may worry about this also injecting grid scale. There are some differences which remain even in the limit phi,A|| = 0 – the inhomogeneous solve includes the solution from the previous/current time step unlike the homogeneous one and we integrate with a starting value of 0 rather than 1 as in the homogeneous problem.

Point 6 and 7 suggests that this may be related to an inconsistency in where we evaluate the source terms compared to the other quantities. However, dropping bakdif to 0 actually indicates we evaluate everything at grid centres, so it is not simply an inconsistency but rather a requirement for some decentering/dissipation on the A|| source. (I can note that the final mode structures are still not particularly smooth even with moderate bakdif). The solution may simply be add bakdif dissipation on A|| source. This is a little unsatisfactory as we know that at least the opt_source = F implementation causes a lot of dissipation and forces one to high theta resolution in order to get linear convergence (issue #91). I think we may get more dissipation/poorer convergence as we are not necessarily consistent in how/what we decentre. For example opt_source = F bakdif's A|| and centres J0 and v||, whilst one might argue we should bakdif v|| A|| J0 together. This is repeated in other places where we grid centre equilibrium quantities (v||, wdrift etc. ) and then multiply these by bakdif'd things. This introduces cross terms not present when centering the entire quantity, e.g. consider centering a*b. If we centre the two quantities and then multiply we get (a1+a2)(b1+b2)/4 = (a1b1 + a2b2 + [a2b1 + a1b2])/4 whilst if we multiply and then centre we get (a1b1 + a2b2)/2. These are clearly different. Writing a2 = a1 + dt da/dt and b2 = b1 + dt db/dt with dt =dtheta we find (a1+a2)(b1+b2)/4 = (a1b1 + a2b2 + [b1 (a1 + dt da/dt) + a1 (b1 + dt da/dt)])/4. In the limit of dt->0 the form (a1+a2)(b1+b2)/4=(a1b1 + a2b2 + [b1 (a1 + dt da/dt) + a1 (b1 + dt da/dt)])/4 reduces to (a1b1 + a2b2)/2 showing that one should converge to the other but clearly the rate depends on da/dt and db/dt. Taking a = A|| and b = v|| , for example, we can see that the error may look like dV||/dt or dA||/dt and one may expect dV||/dt to be particularly significant for trapped particles near the bounce point.