Adding bpar to stella, explicit and implicit algorithms

[mrhardman]

This pull request is intended to deliver a fully functional implementation of the explicit and implicit algorithms using all three fields, phi, apar, and bpar. Starting from the development/apar2 branch where apar is included by introducing the distribution function gbar = g + 2 (Zs/Ts) vths vpa J0 apar, and g = h - Zs J0 phi / Ts, we include bpar. Bpar is included by modifying g to g = h - (Zs/Ts)(J0 phi + 4 (J1/bs) mu (Ts/Zs) bpar).

For the explicit algorithm standard flux tube features are supported (tested linear, collisionless simulations).

For the implicit algorithm collisionless features are supported (tested linear, and nonlinear collisionless simulations).

Initial low-resolution nonlinear simulations have been performed with promising results (no blow-up!).

Linear, collisionless, benchmarks have been carried out against GS2 for an ITG and a KBM example.

Not supported

• perpendicular flow shear for simulations with A|| (The implicit implementation of perpendicular flow shear acts on g, not gbar, meaning that there are (potentially, probably) missing flow-shear terms involving A||). B|| terms are included in the evolved g in the implicit algorithm so there are no missing flow shear terms for B||

• radial variation

• full flux surface

• implicit collisions

• driftkinetic_implicit

• neoclassical terms

• in-code comments: many places refer to g as g = < f >, which is no longer true.

• Support the implicit algorithm for flow shear and A||

• Support the implicit algorithm for collisions.

• Provide electromagnetic fluxes in the diagnostics (are separate variables preferred, or should the existing flux be upgraded to the total one? @mabarnes @GeorgiaActon)

• Provide in-code documentation.

• Provide benchmarks and report documentation.

• (Optional extra) Support neoclassical terms for apar and bpar (@mabarnes where are the equations written down?).

• (Optional extra) Support radial variation.

• (Optional extra) Support full flux surface (@GeorgiaActon I see you have another open PR, so perhaps this is a non-trivial job).

[mrhardman]

@mabarnes @GeorgiaActon The last commits add the bpar contribution to the particle and heat fluxes, but the contribution to the perpendicular part of the momentum flux is not included. To support this stella would need to have the Bessel function d J'_1(x) /d x = (1/2) (J0(x) - J2(x)).

[GeorgiaActon]

Notes on numerical instability:

• Including parallel streaming in isolation seems stable, including mirror in isolation seems stable, but including them both without the presence of drifts can be unstable.

• This numerical instability is worse with increasing dt, decreasing k_y, or decreasing qinp.

• It seems that these parameters are all knobs for the relative amplitudes of the streaming and mirror terms compared to the time derivative.

• Turning off the phi contribution in the streaming removes the instability. However, keeping the phi contribution and reducing the amplitude of the g contribution from streaming kills the instability. The reverse is also true; increasing the amplitude of the g contribution to the streaming equation worsens the instability.

[mrhardman]

Notes on numerical instability:

A further note that all of these experiments were done with zed_upwind = 1.0, demonstrating that even maximum first-order upwinding does not guarantee a stable simulation.

[mrhardman]

An additional note:

Examining the collisionless_mtm_implicit.zip example above, using drifts_implicit = .true. makes the example numerically unstable. This might be consistent with the fact that the implicit algorithm cannot accept a timestep size too large. However, this instability is not repaired by setting include_mirror = .false.. This may suggest that the electromagnetic implementation has addition bugs compared to the electrostatic implementation @GeorgiaActon @mabarnes.

Another check for the instability: turning off the drifts with

&time_advance_knobs
xdriftknob = 0.0
ydriftknob = 0.0
wstarknob = 0.0

removes the time step restraint from the automatic CFL condition and appears to give a stable simulation for the nominal time step (the automatic CFL otherwise reduces the timestep size to 1.0e-3 from 1.0e-2 in the input file).

[mrhardman]

@Scottyujia @mabarnes I have just noticed that the perpendicular flow shear terms connected to apar are likely missing, so perpendicular flow shear is likely not implemented properly in the electromagnetic case. I have updated the tracker at the top of this PR.

[mrhardman]

example_nonlinear_phiaparbpar_2.zip
Example of a very low resolution nonlinear electromagnetic simulation in cyclone base case geometry. This case is definitely under-resolved but shows successful execution of the code. Made with

commit ee47b7ce6c0dacfe00cb165cf7582de1df4b9855 (HEAD -> development/apar2plusbpar, origin/development/apar2plusbpar)

and taking ~1.5 CPU hrs to execute on 16 cores (5.60 mins). Perhaps suitable for "long" automatic tests.

Plots are made with

python3  xarray_post_processing/plot_stella.py example_nonlinear_phiaparbpar_2