Numerical instability for very small ky rhoref in electrostatic stella

[mrhardman]

In trying to understand a numerical instability in a high fidelity linear simulation, I have discovered a numerical instability that affects electrostatic stella. In the attached .zip file is a cyclone base case example which goes numerically unstable after a long period of almost constant amplitude in phi. B is forced to be constant with set_bmag_const = .true., with the intent to remove complications from the mirror term. In the example fprim = 0 and tprim = 0, guaranteeing that this instability is numerical rather than physical. The eigenfunction has a high or grid scale k||. Introducing a larger zed_upwind does delay the instability, but even with zed_upwind = 0.1 the growth of the mode can still be seen.

Sharing this here for future reference. May be of interest to @DenSto @mabarnes @d7919 @rd1042 due to similarities to the famous GS2 numerical instability at low ky and dt (https://bitbucket.org/gyrokinetics/gs2/issues/108/linearly-unstable-zonal-modes), although there the issue seems related also to the electromagnetic terms.

cbc.zip

P.S. for avoidance of doubt, the version of stella used was commit 126b599.

[mrhardman]

As might be expected from reading Appendix A of the stella paper (https://www.sciencedirect.com/science/article/pii/S002199911930066X ) this particular numerical instability can be stabilised for zed_upwind = 1.0, see the attached input file. (Note that set_bmag_const = .true. is commented out -- this flag does not affect this result).

cbc_zed_upwind_1.0.zip

[rd1042]

Just to check - the appendix of the stella paper shows a damped numerical mode (negative growth rate). Numerical instability only occurs if an explicit numerical scheme is used, and the timestep isn't small enough to resolve the spurious mode (a CFL condition).

Given that the input file here indicates a fully implicit approach, I think we're looking at a different instability to that described by the stella paper?

[rd1042]

I've been having a fiddle around with the example input file and have noticed:

• The instability gets worse as dt increases (e.g. from dt=0.01 to 0.1 in the input file provided). This causes instability with both set_bmag_const = .true. and set_bmag_const = .false. This behaviour is the opposite to (https://bitbucket.org/gyrokinetics/gs2/issues/108/linearly-unstable-zonal-modes) (and instabilities I've found in em_stella branch), which get worse as dt falls.
• The instability disappears when include_mirror = .false. in the physics_flags namelist. I had expected this to give the same behaviour as the case with set_bmag_const = .true. (if bmag's constant, the mirror term vanishes), so it's interesting that it doesn't.
• The instability occurs for either implicit mirror option (matrix solve vs. semi-lagrange)

This suggests to me that something's going wrong in the mirror term; a next step would be to examine a simulation which advances the mirror term only.

[rd1042]

The instability disappears when include_mirror = .false. in the physics_flags namelist. I had expected this to give the same behaviour as the case with set_bmag_const = .true. (if bmag's constant, the mirror term vanishes), so it's interesting that it doesn't.

It turns out that set_bmag_const = .true. doesn't result in dbdzed being equal to zero; the flag just forces bmag to a constant value after all the other geometric terms are calculated. In other words, the mirror term is still doing things even if set_bmag_const = .true.

From the comments, set_bmag_const was a temporary flag put in place to test something specific, so I wonder if we still need it (given there are other ways of overwriting bmag.) I'll set up a separate issue to discuss this.

[DenSto]

Does this instability occur with ky = 0 in your set-up? Also, what does the very-long-time behaviour with include_mirror = .false. look like? Is it actually stable, or is it unstable with a much reduced growth rate? I've often seen the latter myself in numerical tests with k_y = 0 (see #59).

EDIT: if this is related to #59, maybe that issue needs to be revisited.

[rd1042]

@DenSto I'll take a look, but some additional (possibly relevant) info:

• I've found phi^2 goes large (from phi^2 (t=0) = 1 to phi^2 (t=20) ~ 1E9) when only the mirror term is included in the simulation (setting (x,y,wstar)driftknob = 0 and include_streaming = .false. ) But if this simulation is run long enough, phi^2 starts to fall. This happens for both implicit options (I haven't tested explicitly). However, g doesn't blow up (tested by writing out maxval(abs(g)) each timestep) . I believe the reason for this behaviour is as follows:
• With only the mirror term included, the distribution function g is just being advected around on the vpa grid (and all z, kx, ky, mu gridpoints are decoupled).
• If the vpa grid extended from -infinity to +infinity, phi wouldn't change at all since g is conserved by the mirror operation. However, the vpa grid is finite, which causes g to "fall of the edge" of the simulation.
• This happens faster for electrons than ions, because of the 1/sqrt(mass) dependency of the mirror term. This means we "lose electrons" relatively quickly and phi become large. Eventually, g(ion) gets lost faster than g(electron) so phi starts falling.

So the result of simulating "mirror term only" is a temporary blowup in phi, caused by the finite vpa grid. I don't know if this would cause issues when the other GKE source terms are included, but it's plausible (and could be remedied by increasing vpa_max).

[rd1042]

Does this instability occur with ky = 0 in your set-up?

Using Michael's cbc input file and setting ky.rho_ref = 0. , kx.rho_ref = 0.1, dt = 0.1, I see a numerical instability manifest after about 100 timesteps. This also occurs when I make the zgrid boundary_option = "self-periodic" and set nperiod=1.

[rd1042]

Also, what does the very-long-time behaviour with include_mirror = .false. look like? Is it actually stable, or is it unstable with a much reduced growth rate?

I've tested up to t=1440 (before the simulation timed out) without seeing an instability. How long did you find you had to run before you see a weakly growing instability?

[mabarnes]

Re: Bob's theory about the mirror term and finite vpa grid, do you find that increasing vpa_max alleviates the initial phi blowup? If the mirror term is all that's in the equation, then after a sufficiently long time the pdf should should be zero for both species, right?

[rd1042]

Re: Bob's theory about the mirror term and finite vpa grid, do you find that increasing vpa_max alleviates the initial phi blowup?

• In the "mirror term only" simulations, increasing vpa_max results in the phi blowup happening more slowly (which I think is consistent with my theory; the same problem is happening, but over a longer timescale).
• In the "full" simulations with dt=0.1, I find increasing vpa_max (from 3.0 to 8.0 to 10.0 to 20.0 at fixed nvgrid=96) slows the blowup but the blowup still occurs; the results between vpa_max=10.0 and vpa_max=20.0 are quite similar which suggests diminishing returns as vpa_max increases.
• Strangely, even with include_mirror=.false. (which is stable), phi^2 jumps around on short timescales; it leaps from phi^2=1 to phi^2=29 on the first timestep, and then oscillates around within a decaying envelope. It would be useful to know what's causing the big jump, as this isn't explained by my mirror term theory.

If the mirror term is all that's in the equation, then after a sufficiently long time the pdf should should be zero for both species, right?

I agree (but haven't run simulations long enough to verify this is what happens).

[rd1042]

Running the same simulation in the branch feature/Davies-linear-EM-terms (commit 39822af ) (keeping beta=0), I find that we still have an instability, but this time, it gets worse as dt gets smaller rather than bigger dt. On the EM branch, I find a the simulation is stable for around 2000 timesteps for dt=5.0 before starting to go unstable (stella master quickly blows up with dt=5.0). Given that the same blowup occurs for either "implicit mirror" option in the master branch, but has a different character in the EM branch, I wonder if the issue is in the formulation of the equations rather than a code bug. To test, we could see if the instability still occurs when running explicitly, though we're restricted by the CFL condition on which dt we can test.

The EM stella branch's fully implicit beta=0 simulation is different to the master branch's fully implicit simulation in 3 big ways:

• The EM branch advances h rather than g=h-Ze/T F_0s . This means that (among other things) that the mirror equation requires a response matrix approach
• The equations and response matrices are contain = <phi^{n+1} - phi^n> rather than phi^{n+1} and phi^n directly
• The drift terms are included in the streaming equation (I wouldn't expect this to make a difference, provided the implementations are bug-free)

In might be interesting to compare this with the maxwellian_inside_zed_derivative option (a third formulation of the equation) but if I understand the code correctly, this isn't fully implemented (it also requires a maxwellian term inside the d/dvpa term of the streaming equation, I think). Feel free to correct me if I've got this wrong though!

[rd1042]

@mabarnes can you remind me why stella uses (vpa, mu) as coordinates rather than (e.g.) (energy, mu)? Does this simplify how trapped particles are simulated, or similar?

[mabarnes]

Yes, it’s because in stellarators with multiple wells the velocity space and field line coordinate would get coupled together in a complicated way if using energy rather than vpa

[mabarnes]

in case it's helpful, I've created a new branch (development/delphi_response), in which one has the option to choose between using phi^n+1 and (phi^n+1-phi^n) when solving for g_inh in parallel streaming. There is surprisingly little that needed changing in the code! It is controlled by the input flag 'use_deltaphi_for_response_matrix' in the &knobs namelist. From my initial tests with CBC at ky=0.1, there is no difference with the default case.

[mabarnes]

In case it is helpful, I have started looking at the CBC with ky=0.001 (using the delphi_response branch, but with options set to mimic the master branch). Here is what I've found so far:

1. If I turn on only implicit streaming and mirror, the code appears to be stable for any time step size chosen.
2. If the magnetic drift is added in explicitly, the code is violently unstable unless a very small time step is taken. This should not be the case based on the 'advection' CFL condition; my suspicion and recollection from when I first wrote the stella paper is that this is the term limiting the time step at small ky.
3. if the magnetic drift is added using the current drifts_implicit=T option, then a slow-growing numerical instability persists for large delt. It can be alleviated but not completely eliminated by setting mirror_semi_lagrange=F.

The only other thing of note is vpa_upwind and zed_upwind being reduced to 0.005.