Electromagnetic case, bakdif nonzero gives error

[billdorland]

Summary

In the course of a study of high-beta stability (ideal ballooning and gyrokinetic) Dylan Langone found a problem with using non-zero values of bakdif. In general, one would like to be able to use this knob for various reasons.

Equilibrium information

The equilibrium that revealed this problem is described in the efit eqdsk format. It is the attached file pierre.efit. All of the input files provided refer to this equilibrium data. The report is focused on the radius rhoc = 0.95. Note that we are using bishop=4, to adjust the beta gradient (and in principle, the magnetic shear, though we did not change that in this case) to be different from the actual value at that radius. An ideal ballooning stability diagram is attached (NoseRho95.png) to illustrate this point. The blue point is where the actual values of s_hat and beta_prime are in the true equilibrium at this radius. The red points are the stability boundary for theta_0 = 0 at this radius. The red points were found with ball. The case we are submitting has the beta gradient enhanced to be double the equilibrium value, so that it is well inside the ideal ballooning unstable region. The fprim and tprim and so forth are all set to be consistent with the enhanced pressure gradient.

Expectation

When running this case with gs2, we expect to find an instability with a growth rate in the range of 1 to 1.5 in the units selected in the input files provided.

Basic observation/report

Dylan was using an executable supplied by me to analyze this equilibrium last year. Unfortunately, I do not have the hash for the version of that source. That is my fault (Dorland). It is at least as new as the subversion revision 4239, and possibly much newer.

We observe that with the old executable, we could use non-zero values of bakdif and obtain reasonable results, in line with basic expectations. With a recent version of the code

(hash = 474cacded62da4bd36bed9a8d59cf5dbcfdd4aab)

a non-value of bakdif produces the wrong result. In the specific case presented here, the growth rate is far too small with the recent source, and the eigenfunctions look incorrect. Visually, the errors in the eigenfunctions have a 2 pi periodicity, suggesting a problem with the wfb and/or the changes in the bakdif implementation for the A_parallel terms. Note that we find problems for various values of collisionality, but in the examples we have no collisions, to isolate the problem.

Details

We are submitting three cases, with an input file and output file for each. (run.in, run.out).

Two of the cases are qualitatively correct. These are r95noneq and r95noneqfexpr0.5). These cases have bakdif = 0., and different values of fexpr.

The third case is qualitatively incorrect. This is r95noneqwrong. The growth rate is roughly two orders of magnitude too small.

Background

A poster from the 2019 Sherwood meeting is also attached. This poster shows the kind of analyses that Dylan was able to perform with the old executable. All results in this poster were very carefully checked before the conference. When we updated to the more recent trunk version and moved to Cori (at NERSC) we found that we needed to set bakdif = 0 or we found problems like the one reported above.

[David Dickinson]

I suspect this is closely related to issue #91 and the changes introduced in commit ab0d9949d370f1cb23327d71bb8e8d22ca810f88

The diagnosis of this issue simply suggested the non-zero bakdif cases with Apar resulted in excessive dissipation – at sufficiently high ntheta the results should converge.

I believe the primary motivation for the addition of bakdif on Apar was to suppress a numerical instability in nonlinear EM runs, I expect this stabilisation was a direct result of the high levels of dissipation this seems to introduce. I think the secondary motivation was to treat all potentials consistently in the source term. It’s not immediately clear to me why bakdif seems to have such a strong effect when applied to Apar , my only thought so far is that we might not be decentring some other terms that are used in the Apar part of the source calculation so it introduces a slight inconsistency (e.g. issue #54).

There should be a workaround currently in setting opt_source = .true. in dist_fn_knobs in the input file – this uses a source calculation that excludes bakdif on Apar as original versions of gs2 did.

Should we consider reverting the introduction of bakdif in the Apar sources? Interestingly the old documentation at http://gyrokinetics.sourceforge.net/wiki/index.php?title=GS2_Input_Parameters&section=16#dist_fn_species_knobs does recommend setting bakdif = 0 for electromagnetic runs but it doesn’t give justification.

[William Dorland]

I vote for not reverting. I think progress on bakdif has been made, and I'd like to hold on to that. I was the source of the advice on bakdif for EM runs long ago, but I was focused on nonlinear EM runs. I had not seen a problem (or I don't remember it anyway) with bakdif in linear EM runs. There is something in the new coding that could possibly be updated a bit more?

You may be right about the dissipation angle, but this is quite an extreme example of how a very small decentering led to a huge change in the solution. I do not believe the bakdif knob should be that sensitive.

[David Dickinson]

Yes it does seem like quite an extreme effect, surprisingly so. Whilst the behaviour in issue #91 is only a factor ~2 difference in the growth rate for moderate resolutions it is also quite an extreme effect when one considers the ntheta required for convergence.

I’m certainly happy to not revert this code as well – I was quite hopeful when it came in that this might help with some difficulties I had been having with certain EM runs.

A postdoc here at York (Michail Anastopoulos Tzanis) wrote a toy code trying to mock up a simple version of the gs2 algorithm to explore this a little while ago but didn’t find a similar effect, which makes me wonder if we’re just currently being a little inconsistent somewhere in the main code.

This is clearly an important issue to resolve.

‌

[David Dickinson]

I’ve confirmed that reverting the changes in ab0d9949d370f1cb23327d71bb8e8d22ca810f88 does cause this problem to go away, indicating it is indeed the Apar related bakdif treatment that has triggered this.

[David Dickinson]

Just to provide a quick update – I’ve started looking at a periodic slab with no equilibrium gradients to see if I can isolate the issue in a simple problem. I can see some issues starting to develop in this case with bakdif and/or electromagnetic terms. This is interesting as for this case almost everything in the equilibrium is constant along the field so spatial decentering might be expected to have less of an impact in some places.

[William Dorland]

Interesting. I should have some time to contribute this weekend. I can try my own ideas or I can work off of yours if you want. I suppose I would need input files from you?