long-wavelength numerical instability in unsheared slab with explicit streaming

[rd1042]

I've tried running a low-res simulation in an unsheared slab in stella (periodic BCs), in the absence of any gradients. I'm finding that when I treat parallel streaming implicitly, I get something stable ( |phi^2| fluctuates but slowly falls over time); however, with explicit streaming the solution blows up; this blowup happens more quickly at lower ky.

I think is related to the electrostatic shear alfven wave problem (described in the appendix of the 2019 stella paper, https://doi.org/10.1016/j.jcp.2019.01.025 ); one expects a numerical instability if the d(phi)/dz term is treated in an upwinded way, and from a look at the code (parallel_streaming.f90), advance_parallel_streaming_explicit calculates d(phi)/dz by calling the subroutine third_order_upwind_zed.

It seems advance_parallel_streaming_explicit used to calculate d(phi)/dz in a centered way, but has been replaced with an upwinded method. I was wondering if anyone knew the reason for this?

Example input files to reproduce instability attached.

input_slab_explicit.zip

[DenSto]

Thanks for bringing this up! I seem that I may have made a goof here... my original intention was probably to use the get_dgdz_variable function here, which instead sets the amount of upwinding by a parameter specified in the input file (which is usually ~ 0.02 or something similarly small). I didn't realize I had set the derivatives there to full upwinding, or at least I can't remember why I chose to do it that way.

If you change the call back to dgdz_centered (which has been commented out but is still in the file) does this problem go away? If so, if you change the function call to dgdz_variable and use a small amount of upwinding, does it also go away?

[DenSto]

Actually, get_dgdz was already there, I just took it out of a do loop.

[mabarnes]

For this term we need to use central differences for the phi term but can be agnostic about the differencing of g (either full upwinding or something specified by the user). Is this just a matter of uncommenting code?

[DenSto]

I've tried running a low-res simulation in an unsheared slab in stella (periodic BCs), in the absence of any gradients. I'm finding that when I treat parallel streaming implicitly, I get something stable ( |phi^2| fluctuates but slowly falls over time); however, with explicit streaming the solution blows up; this blowup happens more quickly at lower ky.

I think is related to the electrostatic shear alfven wave problem (described in the appendix of the 2019 stella paper, https://doi.org/10.1016/j.jcp.2019.01.025 ); one expects a numerical instability if the d(phi)/dz term is treated in an upwinded way, and from a look at the code (parallel_streaming.f90), advance_parallel_streaming_explicit calculates d(phi)/dz by calling the subroutine third_order_upwind_zed.

It seems advance_parallel_streaming_explicit used to calculate d(phi)/dz in a centered way, but has been replaced with an upwinded method. I was wondering if anyone knew the reason for this?

Example input files to reproduce instability attached.

input_slab_explicit.zip

Which branch are you running this on, anyway? I'm getting different behaviour on my local one (even the implicit blows up, albeit much slower)

[rd1042]

I'm running on master (I've just pulled the latest version); attached is the input & output of running implicitly (after t=60, |phi|^2 = 9.0154E-07 )

I think it's just a case of uncommenting code, giving it a try now.

input_slab_implicit.zip

[DenSto]

So here's an update: stream_implicit for me still fails and explodes. I've changed the dgdz in parallel_streaming to centered and the numerical instability goes away. I then changed it to variable and the numerical instability was absent for zed_upwind = 0.02, but appeared again around zed_upwind = 0.1. My inclination is to just use dgdz_variable for everything (and to be more concrete, rename get_dgdz_variable -> get_dgdz and remove the other ones).

It'd be nice to figure out why my implicit solve is failing, however.

[rd1042]

Curious; I tried changing dgdz in parallel_streaming to centered (just for the phi term, and then for phi and g), and running the code explicitly I still get a blowup. Just realised that my remote is the "old" stella repo, so I'll try switching to the new one and repeat.

I think we probably want to make it such that d(phi)/dz always centered (to avoid a known numerical instability), but allow dg/dz to be upwinded as desired by the user

[rd1042]

To confirm, I get the same results as reported by me above when working with the latest repo (explicit always unstable, implicit always stable for the test cases posted on this thread). Just to check Denis, are you running the same sims, also on master? Seems like we both get a blowup, but in different situations.

[DenSto]

Hmm... in your zip file you sent there's only a slab.geometry file, while the input file takes unsheard_slab.geometry, and so I tweaked the input file. Should there be an unsheard_slab.geometry file?

[rd1042]

Oh dear; my mistake - there should be an unsheard_slab.geometry. Corrected inputs attached.

input_slab_explicit.zip input_slab_implicit.zip

[DenSto]

OK, that makes more sense.

My implicit scheme is now stable. The explicit scheme is stable if I use centered (I forgot to point out: zonal_mode(iky) in the derivative function call in dgdz_centered needs to be changed to periodic!), but the results don't agree with the implicit method... the saturated value of |phi^2| is an order of magnitude higher in the explicit scheme. If I use a dgdz_variable, then the scheme is only stable if I use a very small value of zed_upwind (something like 0.0001)

The fact that the implicit and explicit schemes don't agree worries me. Even if running with slab.geometry is pathological (because I'm not sure how it is generated), I would hope they would agree there too.

[DenSto]

For reference, purple is explicit, green is implicit.

[rd1042]

Thanks - applying the change in dgdz_centered, I agree with you. It looks like this mode(ky=0.1) is weakly unstable when we run explicitly, and stable when implicit. I'm also a bit worried, since when I change the d/dz of the phi term to centered but not the dg/dz term, the explicit simulation is unstable (as Michael says, one would expect the numerical instability to be agnostic to how dg/dz treated. In addition, while ky=0.1 is weakly unstable with what we've done, dropping ky further (say, 0.01) gives me instability.

For FYI, slab.geometry is a half-cooked geometry file (it's not actually slab, so the drifts are finite, and the geometry isn't periodic), so might be harder to debug its behaviour

[DenSto]

Hmm... I may look at standard CBC and see if I can get this sort of disagreement between implicit and explicit schemes for very small ky.

Even without any upwinding, the explicit scheme can still suffer numerical instability, but it scales must less strongly with 1/k_y. Perhaps for such small values of ky (like 0.01), this is the numerical instability that we're hitting.

[rd1042]

From a look at the phi, it looks like the solution's getting big at the boundary (example plots attached, 2 different timestep sizes); could be an issue with the BCs / how the finite difference routines handle the BCs. I'll try writing some tests for finite_differences, which should help nail down if this is the case.

[DenSto]

Here's a problem I found in the function fill_zed_ghost_zones in extended_zgrid.f90:

subroutine fill_zed_ghost_zones (it, iseg, ie, iky, g, gleft, gright)

use zgrid, only: nzgrid

implicit none

integer, intent (in) :: it, iseg, ie, iky
complex, dimension (:,:,-nzgrid:,:), intent (in) :: g
complex, dimension (:), intent (out) :: gleft, gright

! stream_sign > 0 --> stream speed < 0

if (iseg == 1) then
   if (periodic(iky)) then
      gleft = g(iky,ikxmod(iseg,ie,iky),iz_up(iseg)-2:iz_up(iseg)-1,it)
   else
      gleft = 0.0
   end if
else
   gleft = g(iky,ikxmod(iseg-1,ie,iky),iz_up(iseg-1)-2:iz_up(iseg-1)-1,it_left(it))
end if

if (nsegments(ie,iky) > iseg) then
   ! connect to segment with larger theta-theta0 (on right)
   gright = g(iky,ikxmod(iseg+1,ie,iky),iz_low(iseg+1)+1:iz_low(iseg+1)+2,it_right(it))
else
   ! apply periodic BC where necessary and zero BC otherwise
   if (periodic(iky)) then
      gright = g(iky,ikxmod(iseg,ie,iky),iz_low(iseg)+1:iz_low(iseg)+2,it)
   else
      gright = 0.0
   end if
end if

end subroutine fill_zed_ghost_zones

This assumes the distribution is periodic on any 2 pi theta segment, but actually for ky .ne. 0, we want it to be periodic over the whole domain, right? In that case we want

 subroutine fill_zed_ghost_zones (it, iseg, ie, iky, g, gleft, gright)

use zgrid, only: nzgrid

implicit none

integer, intent (in) :: it, iseg, ie, iky
complex, dimension (:,:,-nzgrid:,:), intent (in) :: g
complex, dimension (:), intent (out) :: gleft, gright

integer :: nseg

! stream_sign > 0 --> stream speed < 0

nseg = nsegments(ie,iky)

if (iseg == 1) then
   if (periodic(iky)) then
      gleft = g(iky,ikxmod(iseg,ie,iky),iz_up(nseg)-2:iz_up(nseg)-1,it)
   else
      gleft = 0.0
   end if
else
   gleft = g(iky,ikxmod(iseg-1,ie,iky),iz_up(iseg-1)-2:iz_up(iseg-1)-1,it_left(it))
end if

if (nsegments(ie,iky) > iseg) then
   ! connect to segment with larger theta-theta0 (on right)
   gright = g(iky,ikxmod(iseg+1,ie,iky),iz_low(iseg+1)+1:iz_low(iseg+1)+2,it_right(it))
else
   ! apply periodic BC where necessary and zero BC otherwise
   if (periodic(iky)) then
      gright = g(iky,ikxmod(iseg,ie,iky),iz_low(1)+1:iz_low(1)+2,it)
   else
      gright = 0.0
   end if
end if

end subroutine fill_zed_ghost_zones

The question then is, do we want zonal modes to remain periodic on a 2 pi segment, or on the entire z domain? If so, this will have to be modified further. Anyway, see if this makes a difference.

[DenSto]

If we want to keep the old method for zonal modes, then we'll need

subroutine fill_zed_ghost_zones (it, iseg, ie, iky, g, gleft, gright)

use zgrid, only: nzgrid
use kt_grids, only: zonal_mode

implicit none

integer, intent (in) :: it, iseg, ie, iky
complex, dimension (:,:,-nzgrid:,:), intent (in) :: g
complex, dimension (:), intent (out) :: gleft, gright

integer :: nseg

! stream_sign > 0 --> stream speed < 0

nseg = nsegments(ie,iky)

if (iseg == 1) then
   if (periodic(iky)) then
     if (zonal_mode(iky)) then
       gleft = g(iky,ikxmod(iseg,ie,iky),iz_up(iseg)-2:iz_up(iseg)-1,it)
     else
       gleft = g(iky,ikxmod(iseg,ie,iky),iz_up(nseg)-2:iz_up(nseg)-1,it)
     endif
   else
      gleft = 0.0
   end if
else
   gleft = g(iky,ikxmod(iseg-1,ie,iky),iz_up(iseg-1)-2:iz_up(iseg-1)-1,it_left(it))
end if

if (nseg > iseg) then
   ! connect to segment with larger theta-theta0 (on right)
   gright = g(iky,ikxmod(iseg+1,ie,iky),iz_low(iseg+1)+1:iz_low(iseg+1)+2,it_right(it))
else
   ! apply periodic BC where necessary and zero BC otherwise
   if (periodic(iky)) then
      if (zonal_mode(iky)) then
        gright = g(iky,ikxmod(iseg,ie,iky),iz_low(iseg)+1:iz_low(iseg)+2,it)
      else
        gright = g(iky,ikxmod(iseg,ie,iky),iz_low(1)+1:iz_low(1)+2,it)
      endif
   else
      gright = 0.0
   end if
end if

end subroutine fill_zed_ghost_zones

[DenSto]

Yeah, with these changes the two simulations now broadly agree,

We need to make a decision on the zonal flows.

[mabarnes]

The periodic BCs that were implemented in the code were only intended to be used for cases with zero magnetic shear. For such cases, each 2pi segment should be self-periodic for all ky. Was the bug for a case with finite magnetic shear and periodic BCs? I suppose periodic vs zero incoming BCs should agree if the domain is extended far enough along the field line in this case. If they don't, I'm not sure that periodic is a sensible option.

As for zonal flow, as ky=0, I think each 2pi segment should be self-periodic. I.e., the change in radial wavenumber due to the twist-and-shift BC should be zero, so there should be no connections.

[DenSto]

In that case, is it just incorrect to run with nperiod that's not equal to 1 when using no shear? Otherwise, what's enforcing periodicity for iseg not equal to 1 or nseg? The boundary between segments is going to be fishy one way or another.

As far as the zonal mode goes, does it have to enforced, i.e. the zonal flow is the constant/periodic piece over the zgrid. Is it overkill to enforce 2pi periodicity on the non-zonal component as well? Again, the old method only enforced this 2pi periodicity on the first and last segments.

[DenSto]

Actually, if 2pi periodicity is enforced for all ky for every segment, then there is no reason to run with nperiod not equal to 1, since all the segments would be decoupled anyway

[mabarnes]

If one buys the statistical periodicity argument, then there should be a physical 2pi periodicity in a tokamak. As such, enforcing periodicity after N 2pi turns is not inconsistent but does not seem to enforce the actual constraint we would want.

If nperiod > 1 and running with 'default' parallel BC, i.e., ballooning, then there should only be 1 segment for all ky. The only time that there is > 1 segment should be for twist-and-shift BCs in 'box' mode. Here, if magnetic shear is zero, then there should be no kx shift at 2pi boundaries, and so each segment should be self-periodic. This should result in again having only 1 segment for all ky. The only occasion where we would have more than one connected segment with periodic BCs is if magnetic shear is non-zero -- and it's unclear to me if periodic BCs make sense here.

[mabarnes]

@DenSto Have you observed a case with zero magnetic shear and more than one segment? If so, that sounds like a bug.

[rd1042]

My test case which we've been using was an unsheared slab (the "shear" was set to finite, but geometric quantities overwritten) with nperiod=2, and so more than 1 segment. If one sets the shear to zero, one still sees instability, in the absence of the changes Denis just made.

I guess in a slab geometry, we wouldn't expect physical 2pi periodicity - but perhaps it's better always set nperiod=1 for these for simplicity.

[DenSto]

If nperiod > 1 and running with 'default' parallel BC, i.e., ballooning, then there should only be 1 segment for all ky.

That's right, but the code actually breaks this one segment into 2pi chunks, and so nsegments is actually more than one:

case default

   neigen = nakx ; neigen_max = nakx

   if (.not. allocated(ikx_shift_end)) then
      allocate (ikx_shift_end(neigen_max))
      allocate (ikx_shift(nakx,naky))
   end if
   ikx_shift = 0 ; ikx_shift_end = 0

   if (.not. allocated(nsegments)) then
      allocate (nsegments(neigen_max,naky))
   end if

   ! this is the number of 2pi poloidal segments in the extended theta domain,
   ! which is needed in initializing the reponse matrix and doing the implicit sweep
   nsegments = 2*(nperiod-1) + 1

and so the bug is that 2pi periodicity was enforced for the first and last subsegments, but not anything in between. This is only a problem for the 'self-periodic' boundary conditinon (not 'default' or 'linked') and I suppose it could be a problem for the ky= 0 if nperiod is larger than 1, but that's typically not done in nonlinear simulations, and I suppose the zonal mode isn't included in linear simulations with extended modes.

[mabarnes]

My fault -- I had forgotten how nsegments was defined for the 'default' BC. I agree that nsegments > 1 for nperiod > 1

[DenSto]

Again, we could enforce 2pi periodicity in each segment, but then for the 'self-periodic' BC, running with nperiod > 1 doesn't really make any sense, right? (both numerically and physically)

[mabarnes]

Can anyone think of a case where having nperiod > 1 with periodic BCs is sensible? I guess it depends on how the slab geometry is set up; i.e., what is the parallel length scale set by?

[mabarnes]

I suppose the question is whether we force nperiod=1 when using periodic BCs or we use Denis's fix + a similar fix for the zonal modes. Would the latter fix ever be needed for anything of use? If someone wants to run a flux tube (linked) nonlinear simulation with more than one poloidal turn (nperiod > 1), then I'm not even sure that the code is setup to handle it. I.e., I don't recall nperiod being used if 'linked' BCs are employed.

[DenSto]

I'm leaning toward setting nperiod=1 if periodic BCs are used. The problem there is that the zgrid is set up before the condition shat = 0 is checked in stella_geometry. I suppose this can be fixed by breaking up the read_parameters and initialisation in the geometry and zgrid modules and calling those first. This is done already for a few other files.

[DenSto]

Arg, that doesn't work because shat is set by the geometry submodules, like miller and vmec. Maybe one can just print out a warning or perhaps abort altogether?

[rd1042]

Can anyone think of a case where having nperiod > 1 with periodic BCs is sensible? I guess it depends on how the slab geometry is set up; i.e., what is the parallel length scale set by?

What about an unsheared stellarator? If ones imagines (say) a W7X case where shear=0, but q is such that we need to make 6 poloidal turns to get back to the same point, is this a case of setting nfield_periods=6 and nperiod=1? Or would one need to change nperiod?

[DenSto]

OK, I've made changes to extended_zgrid.f90 (the first set of changes that I mentioned), which I think are the most flexible ones, and should allow for low-shear simulations on rational q surfaces.