Addition of complete electromagnetic momentum fluxes

[physycola]

Summary of changes:

The current implementation of the momentum flux in CGYRO is missing some electromagnetic contributions in both the $\delta A_\parallel$ and $\delta B_\parallel$ components. The attached document explains which terms have been neglected, and shows how to write them in a form compatible with the CGYRO source code. Specifically, we show that implementing the previously neglected terms only requires modifying the array jxvec_c in cgyro/src/cgyro_init_arrays.F90. The changes also do not involve any changes to the memory requirements at run time, only during initialisation.

Benchmark results:

These changes have been benchmarked against an equivalent set of momentum-flux diagnostics that were concurrently implemented in GENE. We consider a single linear mode at non-zero $\theta_0$ (PX0) and at finite $\beta_e$. The relevant input files for both codes are attached (input_files.zip). Convergence checks in both the spatial and velocity resolution were performed.

In the first figure below, we plot the growthrate and frequency of the mode as a function of time, showing good agreement between the codes (GENE in blue and CGYRO in red). The eigenfunctions (second figure) show similarly good agreement. In the third figure, we plot the "quasilinear" particle and heat fluxes as a function of time; these are, up to normalisations, defined in the same way as those output directly from a linear CGYRO simulation. In this figure, the solid and dash-dot lines are the electrostatic ($\delta \phi$) and electromagnetic (sum of $\delta A_\parallel$ and $\delta B_\parallel$) components of the fluxes, respectively. The reason that we do not show the $\delta A_\parallel$ and $\delta B_\parallel$ components separately in this case is because GENE does not separate them for the particle and heat fluxes. Finally, the result of the benchmark is shown in the fourth figure, where we plot the contributions to the momentum flux for each field. There is excellent agreement across all channels, where all errors are smaller than the relative (and indeed absolute) errors in the particle and heat fluxes.

Effect of changes:

To demonstrate why these changes are necessary within CGYRO, we run the same linear simulation on commits 740bb2dc811bc1ad38be65ea4b87330995931305 (before the changes) and eb589caf0ced4ead2ccda5f1714eb1561fb7e543 (after the changes). The results of this are shown in the final figure below, where the columns are the contributions from each species (electron, ion, total), while the rows are the contributions from each field ($\delta \phi$, $\delta A_\parallel$, $\delta B_\parallel$). The green lines correspond to a run before the changes, while the red lines correspond to a run after.

As expected, the $\delta \phi$ component of the momentum flux is unchanged, and is dominated by the ion contribution. However, the $\delta A_\parallel$ component is changed quite substantially for the electrons, reversing in sign and almost doubling in magnitude. This is due to the fact that the diagnostics were missing the contribution from the Maxwell stress, which can be significant for the electrons as they carry most of the current. The overall effect of this is an almost doubling of the total momentum flux in the $\delta A_\parallel$ channel. Finally, the $\delta B_\parallel$ channel sees a significant reduction in overall momentum flux, coming from the suppression of the electrons. This is because without the added terms in this channel, the $\delta B_\parallel$ momentum flux is overestimated by a factor of $(k_\perp \rho_s)^{-2}$, which is significant for electrons at the wavenumber considered ($k_y \rho_D = 0.3$). An explanation for this is provided in the attached document.

Overall, these changes mean that the momentum flux calculated in CGYRO is now a complete description of the turbulent momentum flux generated by the fluctuations, capturing the physics associated with the Maxwell stress and the transport of momentum along perturbed field lines that was previously absent.

Figures

Growth rate and frequency:

Eigenfunctions:

Particle and heat fluxes:

Momentum flux:

Momentum flux comparison:

[physycola]

The regression tests pass on the first commit in the pull request, but not on the second. The only regression tests that fail are those with finite BETAE_UNIT (reg03, reg06, reg08, reg09, reg15, reg18, reg19, reg20), and the differences are in the fourth significant figure due to how low BETAE_UNIT is in these cases. This is somewhat to be expected, as it is in these cases that the output of the fluxes will be changed. Running the tests with reset_flag=1 yields passes on all of them.

[physycola]

These changes have been implemented separately in commit 2567ac9 due to the difficulty of merging branches. This PR will be closed.