Reading and plotting linear outputs#
The following discusses how one is able to obtain and visualise data from linear gyrokinetic
simulations within pyrokinetics.
We first need to load in our data from the desired directory using a Pyro object.
Using GS2 as our example code:
from pyrokinetics import Pyro, template_dir
# Point to GS2 input file
gs2_template = template_dir / "outputs/GS2_linear/gs2.in"
# Load in file
pyro = Pyro(gk_file=gs2_template, gk_code="GS2")
# Load in GS2 output data
pyro.load_gk_output()
data = pyro.gk_output
Saving and loading NetCDF files#
For large simulations, reading raw gyrokinetic output files can be slow and memory-intensive,
particularly when additional post-processing steps are required. To improve performance,
GKOutput objects can be written to a NetCDF file and subsequently reloaded directly.
This workflow can significantly reduce loading time when repeatedly analysing the same simulation, and can also help when working under memory constraints (for example on laptops or shared HPC login nodes).
from pyrokinetics import Pyro, template_dir
path = template_dir / "outputs" / "CGYRO_linear"
# Load from native gyrokinetic output
pyro = Pyro(gk_file=path / "input.cgyro")
pyro.load_gk_output()
# Save processed dataset to NetCDF
pyro.gk_output.to_netcdf("linear_cgyro.nc")
# Later: reload directly from NetCDF (much faster)
new_pyro = Pyro(gk_file=path / "input.cgyro")
new_pyro.load_gk_output(netcdf_file="linear_cgyro.nc")
When a NetCDF file is provided via the netcdf_file argument, the expensive parsing of
raw simulation output is skipped and the dataset is reconstructed directly from disk.
This approach is strongly recommended when:
Working with large parameter scans
Performing repeated post-processing or plotting
Operating under limited memory availability
Sharing processed datasets with collaborators
Here we have loaded up the gk_output object, which is an xarray Dataset. More information about this can be obtained by printing:
>>> print(data)
<pyrokinetics.GKOutput>
(Wraps <xarray.Dataset>)
Dimensions: (theta: 73, kx: 1, ky: 1, time: 21, field: 1, species: 2,
energy: 8, pitch: 37, flux: 3)
Coordinates:
* time (time) float64 0.0 1.01 2.02 3.03 ... 17.17 18.18 19.19 20.2
* kx (kx) float64 0.0
* ky (ky) float64 0.07
* theta (theta) float64 -9.425 -9.25 -9.076 ... 9.076 9.25 9.425
* energy (energy) float64 0.004047 0.1044 0.5516 ... 4.739 5.936 7.25
* pitch (pitch) float64 0.003676 0.01921 0.04652 ... 1.182 1.194
* species (species) <U8 'ion1' 'electron'
* field (field) <U3 'phi'
* flux (flux) <U8 'particle' 'heat' 'momentum'
Data variables:
phi (theta, kx, ky, time) complex128 [([tref] / [mref])**(0.5) * [mref] / [bref_B0])·tref_electron/e/lref_minor_radius] ...
particle (field, species, ky, time) float64 [([tref] / [mref])**(0.5)·([tref] / [mref])**(0.5) * [mref] / [bref_B0])²·nref_electron/lref_minor_radius²/rad] ...
heat (field, species, ky, time) float64 [([tref] / [mref])**(0.5)·([tref] / [mref])**(0.5) * [mref] / [bref_B0])²·nref_electron·tref_electron/lref_minor_radius²/rad] ...
momentum (field, species, ky, time) float64 [([tref] / [mref])**(0.5) * [mref] / [bref_B0])²·nref_electron·tref_electron/lref_minor_radius/rad] ...
growth_rate (kx, ky, time) float64 [lref_minor_radius/([tref] / [mref])**(0.5)] ...
mode_frequency (kx, ky, time) float64 [lref_minor_radius/([tref] / [mref])**(0.5)] ...
eigenvalues (kx, ky, time) complex128 [lref_minor_radius/([tref] / [mref])**(0.5)] ...
eigenfunctions (field, theta, kx, ky, time) complex128 [] (4.22038619873...
Attributes: (12/14)
linear: True
gk_code: GS2
input_file: &kt_grids_knobs\n grid_option = 'single'\n/\n\...
attribute_units: {}
title: GKOutput
software_name: Pyrokinetics
... ...
object_uuid: 7bb419cd-8e4d-4465-b075-8baf37c18d34
object_created: 2023-08-07 16:47:06.639348
session_uuid: cba1d994-3318-4e40-8c4d-c41e042a2b3e
session_started: 2023-08-07 16:46:59.397799
netcdf4_version: 1.6.1
growth_rate_tolerance: <xarray.DataArray 'time' ()>\n<Quantity(0.0564343...
where we can see the different data variables available, including their dimensionality, coordinates and units.
To read our desired variable, we use the syntax data["Data_variable"]:
# Get eigenvalues
eigenvalues = data["eigenvalues"]
growth_rate = data["growth_rate"]
mode_freq = data["mode_frequency"]
For linear simulations, one tends to only have a single ky and kx, and thus
data variables such as growth_rate and mode_frequency are essentially 1D
functions of time. These can be plotted using plot (see xarray’s Plotting for further details):
mode_freq.plot(x="time")
plt.show()
For data variables with higher dimensions, indexing can be performed using the standard
xarray dataset methods, such as .sel and .isel. For example, to plot the phi
eigenfunction at the final time point as a function of theta:
# Plot eigenfunction
phi_eig = np.real(data["eigenfunctions"].sel(field="phi").isel(time=-1))
phi_eig.plot(x="theta", label="Real")
phi_i_eig = np.imag(data["eigenfunctions"].sel(field="phi").isel(time=-1))
phi_i_eig.plot(x="theta", label="Imag")
plt.legend()
plt.show()
Similarly for the linear fluxes, one can again specify the coordinates for the desired data. For example, to plot the electrostatic ion energy fluxes:
# Plot ion energy flux
ion_flux = data["heat"].sel(field="phi", species="ion1").sum(dim="ky")
ion_flux.plot()
plt.show()
And analogously for the field data, for example looking at
the magnitude of the phi fluctuations at \(\theta = 0.0\):
# Plot phi
phi = data["phi"].sel(theta=0.0, method="nearest").isel(ky=0).isel(kx=0)
phi = np.abs(phi)
phi.plot.line(x="time")
plt.yscale("log")
plt.show()
Details regarding normalisations and units can be found in Normalisation conventions.