Numerical instability of low-kx zonal modes

[Scottyujia]

When I run the attached linear case with kinetic electrons and zero beta, I noticed that the (kx=0.05, ky = 0) mode grows over time.

I have tried the following cases:

Unstable zonal mode:

1. zed_upwind = 0, include_mirror = .true. (semi-lag., implicit or explicit)
2. zed_upwind = 0.02, include_mirror = .false.

Stable zonal mode:
3. zed_upwind = 0, include_mirror = .false.

From 1 and 3, it looks like the problem is with the mirror terms, but from 2 and 3, the problems seems to be the streaming bit. I am confused, maybe I have done something wrong in the input?

I noticed similar issues raised by @DenSto, but I did not seem to find the solution?

cbc_stella.in.zip

[GeorgiaActon]

Hi Yujia,

I had a look at this and I think this problem is related to the initial condition. There are 2 options you have available to remove this problem:

1. Set initialise_distribition_option = 'rh' -> This is an option for initialising the zonal mode only and looking at its oscillation and residual behaviour. In this case, everything behaves as expected.
2. On line 566 you can change abs(sec(is)%z) to just spec(is)%z (i.e. remove the absolute value).

I will look into this further and try to come up with a better, more permanent solution, but this might be useful for now!

[GeorgiaActon]

Here purple is your original input file, green is option 1. (using 'rh' for this initialisation) and blue is option 2 (removing the absolute value around the species charge in the Maxwellian distribution)

[Scottyujia]

Thanks very much for your help, Georgia! Using 'rh' for the initial condition has definitely solved the issue for the case I posted above. I just tried it for a case where the zonal mode grows like crazy and it reduced the growth rate by a lot but the zonal mode still grows slowly.

I would like to try option 2, but I am not sure which line 566 you are referring to?

This is for the new case (see attached). phi2 vs t.

input.in.zip

[Scottyujia]

Oh, I see. You mean the init_g.f90 file.

[GeorgiaActon]

Sorry yes, in the new pull request of Stella (Stella-v1.0) in the file STELLA_CODE/arrays/initialise_distribution_function.f90, line 566

[Scottyujia]

Option 2 does not seem to work for the new case, either:

phi2_vs_t

[Scottyujia]

In case you are thinking about up-windings. I have tried setting vpa and zed up-windings to zero and only use hyper dissipation in vpa and zed. It does not get rid of the slow growth.

[GeorgiaActon]

hmm - I have just tried re-running and I don't see any growth with all default unwinding values. Are there any differences that you can see in our input files (mine is attached below)?

Here is the plot I get with that input file, incorporating the change to line 566 in the stella-v1.0 version. I will run it longer just to check though.

input.in.zip

[Scottyujia]

The input you attached is for CBC. The plots I showed above are for the new one which is different:
input.in.zip

[Scottyujia]

Latest update: The problem was solved by using initial condition `rh' and implicit_drifts = .true.

[GeorgiaActon]

If you run the initial condition 'default' with implicit_drifts = .true. does this also work (with fix 2. above)?

[Scottyujia]

I can confirm that it also works. Let me summarise the findings so far for the benefit of others reading this issue:

Explicit drifts:
No numerical instability for the zonal mode if dt is set to be small enough (on the order of 10^(-3)). Previously I thought the zonal mode was linearly unstable, but as Georgia pointed out, it was to do with starting with the wrong initial condition. The mode needs to evolve for some time before it levels off or decay. If Georgia were to continue the purple curve (the one with the usual default initial condition) on the top figure, it should tend to a particular value and then decrease, similar to this:

It is worth pointing out that increasing zed_upwind will make the zonal mode grow to larger amplitude before it decays in the case of standard default initialisation. I personally would recommend setting zed_upwind to be as small as possible, such as 0.002 (setting it to exactly zero makes certain cases numerically unstable unless zed dissipation D_zed is set to be very large ~1 or 10, which is not ideal), as it can cause other numerical issues in the electromagnetic cases. This is discussed in detail in issue #206.

Implicit drifts:
dt can be chosen large in this case (~10^(-2)) and no forever growing zonal mode has been found across various cases and initial conditions. So far, it has been the most robust and practically useful option to avoid the zonal numerical instability.