Bakdif on Apar leads to very slow convergence in ntheta

[David Dickinson]

With the attached case a very large ntheta is required for the growth rate and frequency to converge (well in excess of 128) which is surprising for a weakly shaped Miller case.

Convergence can be dramatically accelerated by setting bakdif=0.0 for both species and indeed it is often recommended to use bakdif=0.0 for electromagnetic runs. It turns out that the changes in commit ab0d9949d370f1cb23327d71bb8e8d22ca810f88 introduced some bakdif weighting on the apar terms that enter the source calculation and it is exactly these terms that lead to this slow convergence. This can be confirmed by either removing the weighting and recompiling or by running with opt_source = .true. (as it turns out the commit above did not implement these changes in the optimised routine, which is a separate issue worth noting). These weighting terms have been introduced here with a comment noting that it suppresses a numerical instability in electromagnetic non-linear runs. It seems this is likely achieved by (unintentionally) adding a significant numerical dissipation/diffusion to the apar related drives.

‌

[Peter Hill]

Related? #54 #67

‌

Edit: sorry, I see this was already mentioned in #97

[David Dickinson]

Here’s an example showing the growth rate vs ntheta for the standard code (orig) with and without opt_source (T, F) and two curves from opt_source = F with bakdif factors removed on either Apar_m or Apar_p. This suggests the “excess” dissipation is coming from the factor on Apar_m.

‌

‌

[David Dickinson]

The linear source term looks something like

v|| A|| J0 and we form this as something like v||_c A||_c J0_c or v||_c A||_u J0_c where a c subscript indicates a quantity is evaluated at a grid centre and a u subscript indicates the quantity is evaluated at the bakdif weighted location. One may argue that we should instead evaluate (v|| A|| J0)_c or (v|| A|| J0)_u, which are not quite the same. For example

2 (v|| A|| J0)_c = (v||_+ A||_+ J0_+) + (v||_g A||_g J0_g) (where + indicates ig+1 and g indicates ig) whilst

2 (v||_c A||_c J0_c) = (v||_+ + v||_g)/2 * (J0_+ + J0_g)/2 * (A||_+ + A||_g) = (v||_+ A||_+ J0_+)/4 + (v||_g A||_g J0_g)/4 + MIXED TERMS = (v|| A|| J0)_c / 2 + MIXED TERMS so the two forms are not identical. (It is clear that the average of A*B is not equal to the average of A times the average of B). One might perhaps expect this to be a particular issue when any of v||, A|| or J0 vary rapidly in theta. For example v|| may change quickly if bmag changes rapidly.

[David Dickinson]

Replacing v||_c A||_c J0_c (original, opt_source = T) with (v|| A|| J0)_u appears to help reduce the extreme dissipation resulting from bakdif on the Apar terms. The below plot is as above but with four extra curves. “Move V||” is used to note cases in which we move v|| out of the source coefficients and upwind (v|| A||) directly in the source term. Purple shows where we do this but don’t upwind (i.e. calculate at centres, shows a minor difference), brown shows where this is done and we upwind on Apar_p and we very closely recover the case with no upwind on Apar, whilst pink and grey do the same but also upwind apar_m, with grey also including J0 in what is upwinded. These show a level of dissipation between that of the original code with opt_source = F and opt_source = T.

‌