The Gkeyll 2.0 Code: Documentation Home¶
“Magic Chicken Software Framework” – Artificial ‘Intelligence’ view on Gkeyll
“Don’t Panic” – The Hitchhiker’s Guide to the Galaxy
Gkeyll v2.0 (pronounced as in the book “The Strange Case of Dr. Jekyll and Mr. Hyde”) is a computational plasma physics code mostly written in C/C++ and LuaJIT. Gkeyll contains solvers for gyrokinetic equations, VlasovMaxwell equations, and multifluid equations. Gkeyll contains abinitio and novel implementations of a number of algorithms, and perhaps is unique in using a JIT compiled typeless dynamic language for as its main implementation language.
The Gkeyll package contains two major parts: the gkyl simulation framework and the the postgkyl postprocessing package. Here you will find documentation for the full Gkeyll package.
For license and authors see License and Authors.
Installing Gkeyll¶
“Mostly Painless” – Petr Cagas
Installation instructions for both the gkyl simulation framework and the postgkyl postprocessing tool are provided below.
gkyl install¶
There are two options for installing gkyl. One can install directly from the source code, or via the Conda package. Installing directly from the source is the preferred option, as this method gives users more control over the installation process. For many users who will wish to run gkyl on a cluster, which will have clusterbuilt versions of the Message Passing Interface (MPI) for parallel simulations, and potentially other gkyl depedencies, the source build will allow users to set the appropriate paths to the cluster installations of these dependencies.
Installing from source (preferred)¶
To install gkyl from source, first clone the GitHub repository using:
git clone https://github.com/ammarhakim/gkyl
Navigate into the gkyl
directory to begin.
In many cases, an installation of gkyl will involve building most of gkyl’s dependencies only require a modern C/C++ compiler and Python 3. The full list of dependencies is:
 C/C++ compiler with C++17 support (But NOT Clang >= 12.0 provided by Xcode 12)
 Python 3 >=3.6
 MPI compiler with MPI3 support (>=openmpi 3.0 or >=mpich 3.0)
 LuaJIT 2.1.0
 ADIOS 1.13.1 (But NOT >=ADIOS 2.0)
 Eigen 3.3.7
 CUDA Toolkit >=10.2 (if building with GPU support, NOT FULLY SUPPORTED)
The following instructions assume that at minimum the user has both a C/C++ compiler with C++17 support and Python 3.
Installing using “machine files” (recommended as part of installation from source)¶
For systems on which gkyl has been built before, the code can be
built in three steps using scripts found in the machines/
directory.
Install dependencies using a
mkdeps
script from themachines/
directory:./machines/mkdeps.[SYSTEM].sh
where [SYSTEM]
should be replaced by the name of the system you are
building on, such as macosx
or eddy
. By default, installations
will be made in ~/gkylsoft/
.
Even on systems which have installations of gkyl dependencies such as MPI, the mkdeps script must be run first to build other gkyl dependencies such as LuaJIT.
Configure
waf
using aconfigure
script from themachines/
directory:./machines/configure.[SYSTEM].sh
NOTE: Steps 1 and 2 should only need to be done on the first
build, unless one wishes to change the dependencies.
A successful waf
configure, on a system without GPU support, will look like:
bash$ ./machines/configure.macosx.sh
./waf CC=clang CXX=clang++ MPICC=/Users/junoravin/gkylsoft/openmpi/bin/mpicc MPICXX=/Users/junoravin/gkylsoft/openmpi/bin/mpicxx out=build prefix=/Users/junoravin/gkylsoft/gkyl cxxflags=O3,std=c++17 luajitincdir=/Users/junoravin/gkylsoft/luajit/include/luajit2.1 luajitlibdir=/Users/junoravin/gkylsoft/luajit/lib luajitsharedir=/Users/junoravin/gkylsoft/luajit/share/luajit2.1.0beta3 enablempi mpiincdir=/Users/junoravin/gkylsoft/openmpi/include mpilibdir=/Users/junoravin/gkylsoft/openmpi/lib mpilinklibs=mpi enableadios adiosincdir=/Users/junoravin/gkylsoft/adios/include adioslibdir=/Users/junoravin/gkylsoft/adios/lib configure
Setting top to : /Users/junoravin/gkyl
Setting out to : /Users/junoravin/gkyl/build
Checking for 'clang' (C compiler) : clang
Checking for 'clang++' (C++ compiler) : clang++
Setting dependency path: : /Users/junoravin/gkylsoft
Setting prefix: : /Users/junoravin/gkylsoft/gkyl
Checking for LUAJIT : Found LuaJIT
Checking for MPI : Found MPI
Checking for ADIOS : Found ADIOS
Checking for EIGEN : Found EIGEN
Checking for Sqlite3 : Using Sqlite3
Checking for NVCC compiler : Not found NVCC
'configure' finished successfully (0.843s)
Build the code using:
./waf build install
The final result will be a gkyl
executable located in the
~/gkylsoft/gkyl/bin/
directory. Feel free to add this directory
to your PATH
environment variable or create an alias so you can
simply call gkyl
.
Note that if MPI was built as well as the part of the installation,
~/gkylsoft/openmpi/bin/
needs to be added to the PATH
as
well. Finally, on some distributions, it is required to add
~/gkylsoft/openmpi/lib/
to the LD_LIBRARY_PATH
environmental
variable.
Machine files for nonnative systems¶
For systems that do not already have corresponding files in the
machines/
directory, we encourage you to add files for your
machine. Instructions can be found in machines/README.md
.
Installing from source manually¶
If machine files are not available, the dependencies, configuration, and build can be done manually.
The first step is to build the dependencies. Depending on your system, building dependencies can be complicated. On a Mac or Linux machine you can simply run the mkdeps.sh script in the installdeps directory. To build dependencies cd to:
cd gkyl/installdeps
First, please check details by running:
./mkdeps.sh h
On most supercomputers you will likely need to use the system recommended compilers and MPI libraries. In this case, you should pass the appropriate compilers to mkdeps.sh as follows:
./mkdeps.sh CC=cc CXX=cxx MPICC=mpicc MPICXX=mpicxx ....
You should only build libraries not provided by the system. In practice, this likely means LuaJIT, ADIOS and, perhaps Eigen. (Many supercomputer centers at DOE already offer ADIOS builds and should be preferred instead of your own builds). A typical command will be:
./mkdeps.sh buildadios=yes buildluajit=yes buildeigen=yes
(in addition, you may need to specify compilers also).
By default, the mkdeps.sh script will install dependencies in $HOME/gkylsoft directory. If you install it elsewhere, you will need to modify the instructions below accordingly.
Once you have all dependencies installed, you can build gkyl itself by cding to the topdirectory in the source. gkyl uses the Waf build system. You do NOT need to install waf as it is included with the distribution. However, waf depends on Python (included on most systems). Waf takes a number of options. To get a list do
./waf help
There are two build scenarios: first, all dependencies are installed in $HOME/gkylsoft directory, and second, you are using some libraries supplied by your system.
If you have installed all dependencies in the gkylsoft directory you can simply run:
./waf configure CC=mpicc CXX=mpicxx
where CC and CXX are names of the MPI compilers to use. Note that in some cases the full path to the compiler may need to be specified. If the compilers are already in your path, then you can omit all flags.
If you need to use system supplied builds, you need to specify more complex set of paths. Although you can do this by passing options to the waf build script, it is easiest to follow these steps:
 Copy the configurepar.shin script to configurepar.sh
 Modify this script to reflect the locations of various libraries on your machine. In particular, if you are using prebuilt libraries you will likely need to change the information about MPI and ADIOS.
 Run the configurepar.sh script
Once the configuration is complete, run the following command to build and install (note: if you are working on a cluster and using environment modules, you may need to load them at this point):
./waf build install
The builds are created in the ‘build’ directory and the executable is installed in $HOME/gkylsoft/gkyl/bin, unless you specified a different install prefix. The executable can only be run from the install directory [1].
If you need to clean up a build do:
./waf clean
If you need to uninstall do:
./waf uninstall
LuaJIT builds easily on most machines with standard GCC compiler. Often, you may run into problems on older gcc as they do not include the log2 and exp2 functions unless c99 standard is enabled. To get around this, modify the src/Makefile in LuaJIT. To do this, change the line:
CC= $(DEFAULT_CC)
to:
CC= $(DEFAULT_CC) std=gnu99
To build on Mac OS X Mojave and beyond you must set the following env flag:
export MACOSX_DEPLOYMENT_TARGET=10.YY
where YY
is the version number of the OSX operating system.
For example, to build on OS X Mojave the env flag is:
export MACOSX_DEPLOYMENT_TARGET=10.14
and to build on OS X Catalina the env flag is:
export MACOSX_DEPLOYMENT_TARGET=10.15
Installing with Conda¶
The gkyl package is also available to be installed via Conda, although
this gives less flexibility for keeping the code uptodate as gkyl development continues.
Once Conda is installed and added
to the PATH
, gkyl can be obtained with:
conda install c gkyl gkeyll
Note, that this will also install all dependencies into the Conda install directory. Often this may lead to some conflicts, particularly for the MPI installation, specially if there is another version of MPI already located in the system. gkyl should be run using the MPI provided by Conda.
In general, having Conda and sourcebuilt gkyl on the same machine can cause confusion. In that case please use explicit paths to the mpiexec and gkyl executable you wish to use when running simulations.
Troubleshooting¶
Having trouble building? We will try to compile a list of suggestions and common error messages in this troubleshooting site.
Footnotes
[1]  The reason for this is that gkyl is in reality a LuaJIT compiler extended with MPI. Hence, for the compiler to find Lua modules (i.e. gkyl specific code) certain paths need to be set which is done relative to the install location. 
Postgkyl install¶
There are two options for getting Postgkyl. The first option is to get the Conda package and the second is to clone the source code repository. Installing via Conda is preferred for the majority of users as it requires literally a single Conda command. What is more, Conda is already the suggested way of installing all the dependencies. On the other hand, the later option has an advantage of always having the most uptodate version and is generally required for users that want to contribute to the code.
Installing with Conda (preferred for nondevelopers)¶
Postgkyl can be installed with Conda with literally a single command:
conda install c gkyl postgkyl
Note that the flag for the Gkeyll channel, c gkyl
, is required
even for updating. However, it can be permanently added.
conda config add channels gkyl
conda install postgkyl
Updates can be downloaded with:
conda update c gkyl postgkyl
Note
To install a new package, users need the write permission
for the Anaconda directory. If this is not the case (e.g. on a computing cluster), one can either
create a Conda environment (see tip below) or
install Conda into the $HOME
directory.
Tip
To create a Conda environment for postgkyl, use
conda create n pgkylenv python=3
Then activate the environment with
conda activate pgkylenv
and install postgkyl using the commands above (or the ones below to install from source).
After install, one must always have the pgkylenv
environment activated in order to use postgkyl.
Installing from source (preferred for developers)¶
Postgkyl source code is hosted in a GitHub repository. To get Postgkyl running, one first needs to clone the repository and install dependencies.
First, clone the repository using:
git clone https://github.com/ammarhakim/postgkyl
Postgkyl has these dependencies, which are readily available thru Conda:
 Click
 Matplotlib >= 3.0
 NumPy >=1.13
 PyTables
 SciPy
 SymPy
 Bokeh
 Adios
 MessagePack for python (msgpackpython)
All these dependencies can be easily obtained from the Gkeyll Conda channel, via
conda install c gkyl postgkyl onlydeps
Once the dependencies are installed, postgkyl can be installed by navigating into
the postgkyl
repository and running
python setup.py install
python setup.py develop
Note that these commands only ever need to be run once (even if one is modifying source code). Changes to the source code will be automatically included because we have installed in development mode.
Building adiospy from source
Adios can also be built manually from the source code. Note that for the manual
build, Adios needs to be already installed and its bin
directory
added to the PATH
(the default Gkeyll location is
~/gkylsoft/adios/bin/
). The standard location for the
wrapper in the Gkeyll installation is
gkyl/installdeps/adiosx.x.x/wrappers/numpy/
. After navigating to that directory,
build and install adiospy via
make python
python setup.py install
This currently does not work out of the box with the clang
compiler because of a deprecated library. This can be overcome
removing the lrt
flag from the line 33 of the Makefile
. The
edited lines 32 and 33 should look like this:
adios.so:
python setup.py build_ext
This will allow to complete the adiospy build successfully and it has no know consequences for Postgkyl.
Switching from Conda version to repository¶
While the Conda build of Postgkyl is the suggested version for the
majority of users, the source code repository is required for any code
contributions. We should stress that when switching between the
different version, it is strongly advised to remove the other
version. Having both may lead to an unforeseen behavior based on the
relative order of components in the PATH
and PYTHONPATH
.
The Conda version can be uninstalled with:
conda uninstall postgkyl
Note on the Windows Linux Subsystem (WSL)¶
Historically, Gkeyll was supported only for Unixlike operating systems. However, since the Windows 10 Anniversary Update in August 2016, it became possible to run Gkeyll even on Windows although with reduced performance. This changed with the introduction of WSL 2 in Spring/Summer 2020. With WSL 2, the performance on Unix and Windows is indistinguishable.
Requirements:
 For x64 systems: Version 1903 or higher, with Build 18362 or higher
 For ARM64 systems: Version 2004 or higher, with Build 19041 or higher
 For WSL 2 (recommended): Build 18362 and a CPU supporting HyperV virtualization
Enabling and installing WSL¶
First, WSL needs to be enabled. This can be done either using GUI or PowerShell. For the former, go to Turn Windows features on or off in the Control Panel and check Windows Subsystem for Linux at the very bottom.
Alternatively, this can be done in a terminal (running the PowerShell as an Administrator):
dism.exe /online /enablefeature /featurename:MicrosoftWindowsSubsystemLinux /all /norestart
The Linux distribution itself is available thru Microsoft Store. Currently, the following flavors are available:
 Ubuntu 16.04 LTS
 Ubuntu 18.04 LTS
 Ubuntu 20.04 LTS
 openSUSE Leap 15.1
 SUSE Linux Enterprise Server 12 SP5
 SUSE Linux Enterprise Server 15 SP1
 Kali Linux
 Debian GNU/Linux
 Fedora Remix for WSL
 Pengwin
 Pengwin Enterprise
 Alpine WSL
It is strongly recommended to then enable WSL 2 for much better performance. First, enable virtualization in Administrator PowerShell
dism.exe /online /enablefeature /featurename:VirtualMachinePlatform /all /norestart
Download and install the kernel update package. And finally, set WSL 2 as the default (again using PowerShell as Administrator).
wsl setversion <distribution name> 2
wsl setdefaultversion 2
Interaction between Windows and Linux; GUI¶
As of September 2020, WSL2 does not directly support GUI; however, this is currently supposed to be in the works. Until the official release, this can be overcome with a 3rd party Xserver like VcXsrv (this is our recommended option as the other option does not seem to work on some configurations) or Xming. Note that when using VcXsrv, the Disable access control checkbox needs to be marked when setting XLaunch. Otherwise, the X11 forwarding would not work properly.
Finally, the $DISPLAY
environmental variable needs to be set up on
the Linux side. In WSL 1 this was as simple as
export DISPLAY=:0
WSL 2, however, uses virtualization and the forwarding address changes
which each launch of the system. The address is stored in
/etc/resolv.conf. The $DISPLAY
can be set up with:
export DISPLAY=$(cat /etc/resolv.conf  grep nameserver  awk '{print $2}'):0
Another issue might be caused by tools that want to open a browser,
for example, the commonly used jupyterlab
. This can be easily
overcome by launching the server without a browser using jupyterlab
nobrowser
and then copying the address including the token to a
Windows browser.
Finally, it is often useful to access the Linux files from
Windows. Using WLS 2, the Linux root, /
, is located at
\\wsl$\<DISTRIBUTION_NAME>
.
Known issues¶
There is currently a known issue where Windows and Linux clocks might
get desynchronized when the computer sleeps. This might cause issues
with Git and update installation using sudo apt update
. There is a
workaround that works until this issue gets patched and that is
manually calling sudo hwclock s
to manually synchronize the time.
Quickstart¶
Once you have a working installation of Gkeyll, you can get to know the code and learn how to use it by following these quickstart examples.
For additional examples one can browse the regression tests folder in the gkyl code or
repository (gkyl/Regression/
). We also began compiling input files for simulations
reported by publications in this repository.
Lastly, more examples and in depth information about the gkyl Apps can be found in the
gkyl Reference.
A first Gkeyll simulation¶
Gkeyll supports several models and numerous capabilities. Before diving into those, here’s a short overview of an input file, how to run it and plot its results.
Contents
Physics model¶
Consider the collisionless Landau damping of an electrostatic wave in a plasma. We simulate it with the VlasovMaxwell model:
More information about the model can be found in the VlasovMaxwell documentation. In this case we consider an electrostatic wave, and there are subtleties regarding electrostatic simulation using Maxwell’s induction equations. So we initialize the electron density as having a small sinusoidal perturbation:
and initialize the electric field in a manner consistent with the Poisson equation (using quasineutrality of the equilibrium electron and ion densities, \(n_i(x,t=0)=n_0\)):
Input file¶
This simulation is setup using Normalized units for the VlasovMaxwell system in a short Lua input file, which begins with the Preamble:
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()
permitt = 1.0  Permittivity of free space.
permeab = 1.0  Permeability of free space.
eV = 1.0  Elementary charge, or JouleeV conversion factor.
elcMass = 1.0  Electron mass.
ionMass = 1.0  Ion mass.
nElc = 1.0  Electron number density.
nIon = nElc  Ion number density.
Te = 1.0  Electron temperature.
Ti = Te  Ion temperature.
vtElc = math.sqrt(eV*Te/elcMass)  Electron thermal speed.
vtIon = math.sqrt(eV*Ti/ionMass)  Ion thermal speed.
wpe = math.sqrt((eV^2)*nElc/(permitt*elcMass))  Plasma frequency.
lambdaD = vtElc/wpe  Debye length.
 Amplitude and wavenumber of sinusoidal perturbation.
pertA = 1.0e3
pertK = .750/lambdaD
 Maxwellian in (x,vx)space, given the density (denU), bulk flow
 velocity (flowU), mass and temperature (temp).
local function maxwellian1D(x, vx, den, flowU, mass, temp)
local v2 = (vx  flowU)^2
local vtSq = temp/mass
return (den/math.sqrt(2*math.pi*vtSq))*math.exp(v2/(2*vtSq))
end
The Preamble typically consists of a line loading the Gkeyll plasma ‘App’ to be used
(in this case VlasovMaxwell), and a specification of a number of user input parameters
and simple derived quantities. One can also create userdefined functions, like
maxwellian1D
in this case, which may be used in the preamble or in the Gkeyll
App that follows. For this specific simulation the Gkeyll app is created by the following:
plasmaApp = Plasma.App {
tEnd = 20.0/wpe,  End time.
nFrame = 20,  Number of output frames.
lower = {math.pi/pertK},  Lower boundary of configuration space.
upper = { math.pi/pertK},  Upper boundary of configuration space.
cells = {64},  Configuration space cells.
polyOrder = 1,  Polynomial order.
periodicDirs = {1},  Periodic directions.
elc = Plasma.Species {
charge = eV, mass = elcMass,
lower = {6.0*vtElc},  Velocity space lower boundary.
upper = { 6.0*vtElc},  Velocity space upper boundary.
cells = {64},  Number of cells in velocity space.
init = function (t, xn)  Initial conditions.
local x, v = xn[1], xn[2]
return (1+pertA*math.cos(pertK*x))*maxwellian1D(x, v, nElc, 0.0, elcMass, Te)
end,
evolve = true,  Evolve species?
},
ion = Plasma.Species {
charge = eV, mass = ionMass,
lower = {6.0*vtIon},  Velocity space lower boundary.
upper = { 6.0*vtIon},  Velocity space upper boundary.
cells = {64},  Number of cells in velocity space.
init = function (t, xn)  Initial conditions.
local x, v = xn[1], xn[2]
return maxwellian1D(x, v, nIon, 0.0, ionMass, Ti)
end,
evolve = true,  Evolve species?
},
field = Plasma.Field {
epsilon0 = permitt, mu0 = permeab,
init = function (t, xn)  Initial conditions.
local Ex, Ey, Ez = pertA*math.sin(pertK*xn[1])/pertK, 0.0, 0.0
local Bx, By, Bz = 0.0, 0.0, 0.0
return Ex, Ey, Ez, Bx, By, Bz
end,
evolve = true,  Evolve field?
},
}
The Gkeyll App typically consists of three sections:
 Common: a declaration of parameters that control the (configuration space)
discretization, and time advancement. This first block of code in
Plasma.App
may specify the periodic directions, the MPI decomposition, and the frequency with which to output certain diagnostics.  Species: Definition of the species to be considered in the simulation. Each species
gets its own Lua table, in which one provides the velocityspace domain and
discretization of that species (for kinetic models), initial condition, diagnostics,
boundary conditions, and whether to evolve it or not (
evolve
).  Fields: A field table, which tells the App whether to evolve the electric and/or magnetic fields according to the field equations of the model. In this table we also specify the initial condition of the fields.
In some applications other sections of the Plasma.App may be necessary, for example, to specify the geometry.
Finally, an input file concludes with an invocation of the App’s run method:
plasmaApp:run()
Running your first simulation¶
Now that we have a Gkeyll input file (named vmdamp.lua
),
simply run the simulation by typing
gkyl vmdamp.lua
You should see the program printing to screen like this:
wsName:gkyldir gabriel$ gkyl vmdamp.lua
Tue Sep 15 2020 16:16:44.000000000
Gkyl built with b0b8203670c7+
Gkyl built on Sep 14 2020 16:29:40
Initializing PlasmaOnCartGrid simulation ...
** WARNING: timeStepper not specified, assuming rk3
Using CFL number 0.333333
Initializing completed in 0.0629927 sec
Starting main loop of PlasmaOnCartGrid simulation ...
Step 0 at time 0. Time step 0.00727108. Completed 0%
0123456789 Step 276 at time 2.00698. Time step 0.00727174. Completed 10%
0123456789 Step 551 at time 4.00677. Time step 0.00727214. Completed 20%
0123456
Gkeyll prints a number every 1% of the simulation, and a longer message with the total number of time steps taken, the simulation time and the latest time step size every 10% of the simulation. This particular simulation ran in 74 seconds on a 2015 MacBookPro. As it progressed it wrote out diagnostic files.
Plotting¶
In this case we did not request additional diagnostics, so the only ones provided are default ones:
 Distribution functions:
vmdamp_elc_#.bp
andvmdamp_ion_#.bp
.  Electromagnetic fields:
vmdamp_field_#.bp
.  Field energy:
vmdamp_fieldEnergy.bp
.
Fields that are larger (in memory) like the distribution function, get written out
periodically, not every time step. These snapshots (frames) are labeled by the number
#
at the end of the file name.
In order to plot the initial distribution function of the electrons we will use the
Gkeyll postprocessing tool (postgkyl), invoked by the pgkyl
command as follows
pgkyl vmdamp_elc_0.bp interpolate plot
This produces the 2D plot of the initial Maxwellian distribution given below.
We can also examine the electrostatic energy in the simulation. This most clearly
exhibits the wave energy decaying as the collisionless damping takes effect. For
this purpose we use the following postgkyl command (we select
the
xcomponent, and pg_cmdplot can use a log scale, as well as add labels):
pgkyl vmdamp_fieldEnergy.bp select c0 plot logy x 'time' y '$E_x^2$'
resulting in the following figure of the (normalized) electrostatic energy as a function of time
Additional quickstart examples¶
The above example used a VlasovMaxwell simulation to showcase how to setup, run and postprocess a Gkeyll simulation. In addition to VlasovMaxwell there are also Gyrokinetic and (fluid) Moment models. Each of these have slightly different features and ways of using them. Quick examples for each of these are found below:
Vlasov example¶
In this example, we extend the VlasovMaxwell input file shown in A first Gkeyll simulation to simulate a more general kinetic plasma. Because the VlasovMaxwell system of equations is widely applicable in plasma physics, this example is intended to illustrate some of the functionality of the VlasovMaxwell solver a user may desire for production science. For more extensive documentation on all of the available options for VlasovMaxwell simulations, we refer the reader to VlasovMaxwell App: VlasovMaxwell equations on a Cartesian grid.
This simulation is based on the studies of [Skoutnev2019] and [Juno2020] and concerns the evolution of instabilities driven by counterstreaming beams of plasma. This example demonstrates the flexibility of the VlasovMaxwell solver by showing how one extends VlasovMaxwell simulations to higher dimensionality, in this case 2x2v. The input file for this example is also a standard performance benchmark for Gkeyll and timings for this input file with varying grid resolution can be found in the gkyl repo in the Benchmarks/ folder.
Contents
Physics model and initial conditions¶
This example solves the VlasovMaxwell system of equations
in two spatial dimensions and two velocity dimensions (2x2v). For this example, the ions are taken to be a stationary, neutralizing background, and therefore do not contribute to the plasma current \(\mathbf{J}\). The electrons are initialized as two equal density, equal temperature, counterpropagating, Maxwellian beams
The electromagnetic fields are initialized as a bath of fluctuations in the electric and magnetic fields in the two configuration space dimensions,
where \(N=16\) and \(\tilde B_{n_x,n_y}\) and \(\tilde \phi_{n_x,n_y}\) are random amplitudes and phases. The electric fields, \(E_x, E_y\) are initialized similarly. Note that the sum goes from \(N\) to \(N\) so as to initialize phases from 0 degrees to 180 degrees.
Input file¶
The full Lua input file (found here) has three components: the App dependencies, the Preamble, and the App. In the App dependencies section, we load the necessary components of Gkeyll to perform a VlasovMaxwell simulation, as well as any additional functionality we require:

 App dependencies

 Load the VlasovMaxwell App
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()
 Pseudorandom number generator from SciLua package for initial conditions.
 Specific rng is the mrg32k3a (multiple recursive generator) of Ecuyer99.
local prng = require "sci.prng"
local rng = prng.mrg32k3a()
The Preamble to set up the initial conditions is:

 Preamble

 Constants
permitt = 1.0  Permittivity of free space
permeab = 1.0  Permeability of free space
lightSpeed = 1.0/math.sqrt(permitt*permeab)  Speed of light
chargeElc = 1.0  Electron charge
massElc = 1.0  Electron mass
 Initial conditions
 1 = Rightgoing beam; 2 = Leftgoing beam.
nElc1 = 0.5
nElc2 = 0.5
ud = 0.1  Drift velocity of beams
uxElc1 = 0.0
uyElc1 = ud
uxElc2 = 0.0
uyElc2 = ud
R = 0.1  Ratio of thermal velocity to drift velocity
TElc1 = massElc*(R*ud)^2
TElc2 = massElc*(R*ud)^2
vthElc1 = math.sqrt(TElc1/massElc)
vthElc2 = math.sqrt(TElc2/massElc)
k0_TS = 6.135907273413176  Wavenumber of fastest growing twostream mode
theta = 90.0/180.0*math.pi  0 deg is pure Weibel, 90 deg is pure twostream
kx_TS = k0_TS*math.cos(theta)
ky_TS = k0_TS*math.sin(theta)
k0_Weibel = 2.31012970008316  Wavenumber of fastest growing Weibel mode
theta = 0.0/180.0*math.pi  0 deg is pure Weibel, 90 deg is pure twostream
kx_Weibel = k0_Weibel*math.cos(theta)
ky_Weibel = k0_Weibel*math.sin(theta)
kx = k0_Weibel
ky = k0_TS/3.0
perturb_n = 1e8
 Perturbing the first 16 wave modes with random amplitudes and phases.
 Note that loop goes from N to N to sweep all possible phases.
N=16
P={}
for i=N,N,1 do
P[i]={}
for j=N,N,1 do
P[i][j]={}
for k=1,6,1 do
P[i][j][k]=rng:sample()
end
end
end
 Domain size and number of cells
Lx = 2*math.pi/kx
Ly = 2*math.pi/ky
Nx = 16
Ny = 16
vLimElc = 3*ud  Maximum velocity in velocity space
NvElc = 16
 Maxwellian in 2x2v
local function maxwellian2D(n, vx, vy, ux, uy, vth)
local v2 = (vx  ux)^2 + (vy  uy)^2
return n/(2*math.pi*vth^2)*math.exp(v2/(2*vth^2))
end
The Preamble defines the constants in the normalization standard outlined in Normalized units for the VlasovMaxwell system and sets the parameters and perturbations to the wave modes of interest for the study. Note that because the dimensionality of the simulation is now 2x2v, the normalization of the Maxwellian has correspondingly changed from the 1x1v Langmuir wave simulation described in A first Gkeyll simulation.
The App can be further subdivided into a number of sections
plasmaApp = Plasma.App {

 Common

...

 Species

...

 Fields

...
}

 Run application

plasmaApp:run()
The Common section of the App defines input parameters which will be utilized by all solvers in the simulation.
For example, the configuration space extents and number of configuration space cells (lower, upper, cells
), as well as what directions, if any, utilize periodic boundary conditions (periodicDirs
), and how to parallelize the simulation (decompCuts,useShared
).

 Common

logToFile = true,
tEnd = 50.0,  End time
nFrame = 1,  Number of output frames
lower = {0.0,0.0},  Lower boundary of configuration space
upper = {Lx,Ly},  Upper boundary of configuration space
cells = {Nx,Ny},  Configuration space cells
basis = "serendipity",  One of "serendipity", "maximalorder", or "tensor"
polyOrder = 2,  Polynomial order
timeStepper = "rk3s4",  One of "rk2", "rk3", or "rk3s4"
 MPI decomposition for configuration space
decompCuts = {1,1},  Cuts in each configuration direction
useShared = true,  If using shared memory
 Boundary conditions for configuration space
periodicDirs = {1,2},  periodic directions (both x and y)
 Integrated moment flag, compute integrated quantities 1000 times in simulation
calcIntQuantEvery = 0.001,
The Species section of the App defines the speciesspecific inputs for the VlasovMaxwell simulation within a Plasma.Species
table.
For example, the velocity space extents and number of velocity space cells (lower, upper, cells
), the function which prescribes the initial condition, and the types of diagnostics.
More discussion of diagnostic capabilities can be found in VlasovMaxwell App: VlasovMaxwell equations on a Cartesian grid.

 Electrons

elc = Plasma.Species {
charge = chargeElc, mass = massElc,
 Velocity space grid
lower = {vLimElc, vLimElc},
upper = {vLimElc, vLimElc},
cells = {NvElc, NvElc},
 Initial conditions
init = function (t, xn)
local x, y, vx, vy = xn[1], xn[2], xn[3], xn[4]
local fv = maxwellian2D(nElc1, vx, vy, uxElc1, uyElc1, vthElc1) +
maxwellian2D(nElc2, vx, vy, uxElc2, uyElc2, vthElc2)
return fv
end,
evolve = true,
diagnostics = {"M0","M1i","M2ij","M3i","intM0","intM1i","intM2Flow","intM2Thermal"},
},
Note that for this particular simulation the ions are a stationary, neutralizing background that does not contribute to the plasma current, so we only require a species table for the electrons.
The Field section if the final section of the App and specifies the input parameters for the field equation, in this case Maxwell’s equation, in the Plasma.Field
table.
For example, similar to the Plasma.Species
table, the Plasma.Field
table contains the initial condition for the electromagnetic field.

 Field solver

field = Plasma.Field {
epsilon0 = permitt, mu0 = permeab,
init = function (t, xn)
local x, y = xn[1], xn[2]
local E_x, E_y, B_z = 0.0, 0.0, 0.0
for i=N,N,1 do
for j=N,N,1 do
if i~=0 or j~=0 then
E_x = E_x + perturb_n*P[i][j][1]*math.sin(i*kx*x+j*ky*y+2*math.pi*P[i][j][2])
E_y = E_y + perturb_n*P[i][j][3]*math.sin(i*kx*x+j*ky*y+2*math.pi*P[i][j][4])
B_z = B_z + perturb_n*P[i][j][5]*math.sin(i*kx*x+j*ky*y+2*math.pi*P[i][j][6])
end
end
end
return E_x, E_y, 0.0, 0.0, 0.0, B_z
end,
evolve = true,
},
Running the simulation¶
The input file vmtsw2x2v.lua
can be run using the gkyl executable
~/gkylsoft/gkyl/bin/gkyl vmtsw2x2v.lua
assuming gkyl
has been installed in the user’s home directory.
When running this simulation, a user should see the following output
Wed Sep 16 2020 11:38:54.000000000
Gkyl built with a4430cbb5d93
Gkyl built on Sep 16 2020 01:25:31
Initializing VlasovMaxwell simulation ...
Initializing completed in 1.39731 sec
Starting main loop of VlasovMaxwell simulation ...
Step 0 at time 0. Time step 0.0360652. Completed 0%
0123456789 Step 139 at time 5.01307. Time step 0.0360652. Completed 10%
01234
The full screen output can be found here, which includes performance details for the simulation. This example was run with a single core of a 10th gen Intel i9 (Comet Lake) processor. Increasing the resolution to \(32^2 \times 32^2\) and now running the simulation using all 10 cores of the Intel i9 using
~/gkylsoft/openmpi/bin/mpirun n 10 ~/gkylsoft/gkyl/bin/gkyl vmtsw2x2v.lua
we obtain the following performance with useShared=true
and the installed MPI from the Gkeyll build.
Postprocessing¶
The output of this simulation is the following set of files:
 Distribution functions:
vmtsw2x2v_elc_#.bp
.  Electromagnetic fields:
vmtsw2x2v_field_#.bp
.  Velocity moments:
vmtsw2x2v_elc_M0_#.bp
,vmtsw2x2v_elc_M1i_#.bp
,vmtsw2x2v_elc_M2ij_#.bp
, andvmtsw2x2v_elc_M3i_#.bp
.  Field energy:
vmtsw2x2v_fieldEnergy.bp
.  Volume integrated moments:
vmtsw2x2v_elc_intM0.bp
,vmtsw2x2v_elc_intM1i.bp
,vmtsw2x2v_elc_intM2Flow.bp
, andvmtsw2x2v_elc_intM2Thermal.bp
.
Snapshots (frames) are labeled by the number #
at the end of the file name, while volume integrated diagnostics that are computed as a timeseries, such as the field energy, are written out as a single file.
Since nFrame=1
in the input file, the only frames that are output are 0
, corresponding to the initial condition, and 1
, corresponding to the end of the simulation.
Since this simulation has two configuration space dimensions, postgkyl creates pcolor plots when run from the command line with pgkyl
.
We can compare the initial condition and final state of the magnetic field, \(B_z\), (of the \(32^2 \times 32^2\) simulation) in two separate figures with the pgkyl
command:
pgkyl vmtsw2x2v_field_0.bp vmtsw2x2v_field_1.bp interp sel comp 5 plot b fixaspect
The default postgkyl colorbar is sequential and useful for visualizing data such as distribution functions, which will vary from 0 (zero phase space density/no particles) to some number (corresponding to a local increase in phase space density).
However, we can see that the colorbar for the magnetic field varies between roughly equal positive and negative numbers, and thus a diverging colormap may yield a more useful representation of the data.
In addition, we can utilize the flexibility of the interpolate
command to interpolate the discontinuous Galerkin data onto an even finer mesh
pgkyl vmtsw2x2v_field_0.bp vmtsw2x2v_field_1.bp interp i 6 sel comp 5 plot b fixaspect diverging xlabel '$x (d_e) $' ylabel '$y (d_e) $'
where we have now added labels with the normalized units in \(x\) and \(y\).
Note that the default interpolation level for polynomial order 2 is 3 (polyOrder
+ 1).
We can likewise visualize diagnostic moments such as the first velocity moment elc_M1i
pgkyl vmtsw2x2v_elc_M1i_1.bp interp i 6 plot fixaspect diverging xlabel '$x (d_e) $' ylabel '$y (d_e) $'
Note that elc_M1i
has two components due to the fact that this simulation has two velocity dimensions, and both components are visualized when this pgkyl
command is utilized.
The left plot is the \(v_x\) velocity moment and the right plot is the \(v_y\) velocity moment.
Further details on the diagnostics available and their definitions can be found in VlasovMaxwell App: VlasovMaxwell equations on a Cartesian grid.
We can also visualize the distribution function from this simulation. However, for this simulation the distribution function if fourdimensional, two configuration space and two velocity space dimensions. Postgkyl offers a number of options for downselecting the data to be more amenable to visualizing. For example, we can readin a subset of the data and visualize the distribution function in velocity space \(v_xv_y\) in the lower left corner of the domain
pgkyl vmtsw2x2v_elc_1.bp z0 0 z1 0 interp i 6 sel z0 0.0 z1 0.0 plot xlabel '$v_x (v_{th_e}) $' ylabel '$v_y (v_{th_e}) $' vmin 0.0
Note that the immediate z0 0 z1 0
tells postgkyl to read in only the first \(x\) and \(y\) configuration space grid cells (while still reading in all of velocity space).
Because we are then interpolating the data onto a finer mesh, the data is still four dimensional so we pass the abbreviated select command sel
to finally downselect to the lowerleft corner of the configuration space domain.
These selective readin commands are vital for very large arrays where the cost in memory and CPU time can be quite large to readin and manipulate the data structure of interest.
Alternatively, if we do want to readin the whole array, we can perform other manipulations to the distribution function such we can still easily visualize the data.
For example, we can use the integrate
command to integrate the distribution function over \(x\) and \(v_x\) to produce a \(yv_y\) plot of the electron distribution function.
pgkyl vmtsw2x2v_elc_1.bp interp integrate 0,2 plot xlabel '$y (d_e) $' ylabel '$v_y (v_{th_e}) $' vmin 0.0
Finally, since we performed this simulation at two different resolutions, and interesting diagnostic to look at is a comparison of integrated quantities between the two simulations.
For ease of plotting we have moved the data from the two simulations to two different folders, res1
(\(16^2 \times 16^2\)) and res2
(\(32^2 \times 32^2\)).
Here, we are being agnostic on what a user might have named these two different simulations and labeling them ourselves with postgkyl.
pgkyl res1/*fieldEnergy.bp l '$16^2 \times 16^2$' res2/*fieldEnergy.bp l '$32^2 \times 32^2$' select comp 5 plot logy xlabel '$t (\omega_{pe}^{1})$' ylabel '$\int B_z^2$' f0
References¶
[Skoutnev2019]  Skoutnev, V., Hakim, A., Juno, J., & TenBarge, J. M. (2019). “TemperatureDependent Saturation of WeibelType Instabilities in Counterstreaming Plasmas”, Astrophysical Journal Letters, 872, (2). https://doi.org/10.3847%2F20418213%2Fab0556 
[Juno2020]  Juno, J., Swisdak, M. M., TenBarge. J. M., Skoutnev, V., & Hakim, A. “Noiseinduced magnetic field saturation in kinetic simulations”, Journal of Plasma Physics, 86, (4). https://doi.org/10.1017/S0022377820000707 
Gyrokinetic example¶
The gyrokinetic system is used to study turbulence in magnetized plasmas. Gkeyll’s Gyrokinetic App is specialized to study the edge/scrapeofflayer region of fusion devices, which requires handling of large fluctuations and openfieldline regions. In this example, we will set up a gyrokinetic problem on open magnetic field lines (e.g. in the tokamak scrapeoff layer). Using specialized boundary conditions along the field line that model the sheath interaction between the plasma and conducting end plates, we will see that we can model the selfconsistent formation of the sheath potential. This simple test case can then be used as a starting point for full nonlinear simulations of the tokamak SOL.
Contents
Input file¶
The full Lua input file (gksheath.lua) for this simulation is a bit longer than the one in A first Gkeyll simulation, but not to worry, we will go through each part of the input file carefully.
To set up a gyrokinetic simulation, we first need to load the
Gyrokinetic
App package and other dependencies. This should be done
at the top of the input file, via

 App dependencies

local Plasma = require("App.PlasmaOnCartGrid").Gyrokinetic()  load the Gyrokinetic App
local Constants = require "Lib.Constants"  load some physical Constants
Here we have also loaded the Constants
library, which
contains various physical constants that we will use later.
The next block is the Preamble, containing input paramters and simple derived quantities:

 Preamble

 Universal constant parameters.
eps0 = Constants.EPSILON0
eV = Constants.ELEMENTARY_CHARGE
qe = eV  electron charge
qi = eV  ion charge
me = Constants.ELECTRON_MASS  electron mass
mi = 2.014*Constants.PROTON_MASS  ion mass (deuterium ions)
 Plasma parameters.
Te0 = 40*eV  reference electron temperature, used to set up electron velocity grid [eV]
Ti0 = 40*eV  reference ion temperature, used to set up ion velocity grid [eV]
n0 = 7e18  reference density [1/m^3]
 Geometry and magnetic field parameters.
B_axis = 0.5  magnetic field strength at magnetic axis [T]
R0 = 0.85  device major radius [m]
a0 = 0.15  device minor radius [m]
R = R0 + a0  major radius of flux tube simulation domain [m]
B0 = B_axis*(R0/R)  magnetic field strength at R [T]
Lpol = 2.4  device poloidal length (e.g. from bottom divertor plate to top) [m]
 Parameters for collisions.
nuFrac = 0.1  use a reduced collision frequency (10% of physical)
 Electron collision freq.
logLambdaElc = 6.6  0.5*math.log(n0/1e20) + 1.5*math.log(Te0/eV)
nuElc = nuFrac*logLambdaElc*eV^4*n0/(6*math.sqrt(2)*math.pi^(3/2)*eps0^2*math.sqrt(me)*(Te0)^(3/2))
 Ion collision freq.
logLambdaIon = 6.6  0.5*math.log(n0/1e20) + 1.5*math.log(Ti0/eV)
nuIon = nuFrac*logLambdaIon*eV^4*n0/(12*math.pi^(3/2)*eps0^2*math.sqrt(mi)*(Ti0)^(3/2))
 Derived parameters
vti = math.sqrt(Ti0/mi)  ion thermal speed
vte = math.sqrt(Te0/me)  electron thermal speed
c_s = math.sqrt(Te0/mi)  ion sound speed
omega_ci = math.abs(qi*B0/mi)  ion gyrofrequency
rho_s = c_s/omega_ci  ion sound gyroradius
 Simulation box size
Lx = 50*rho_s  x = radial direction
Ly = 100*rho_s  y = binormal direction
Lz = 4  z = fieldaligned direction
This simulation also requires a source, which models plasma crossing the separatrix. The next part of the Preamble initializes some source parameters, along with some functions that will be used later to set up the source density and temperature profiles.
 Source parameters
P_SOL = 3.4e6  total SOL power, from experimental heating power [W]
P_src = P_SOL*Ly*Lz/(2*math.pi*R*Lpol)  fraction of total SOL power into flux tube domain [W]
xSource = R  source peak radial location [m]
lambdaSource = 0.005  source radial width [m]
 Source density and temperature profiles.
 Note that source density will be scaled to achieve desired source power.
sourceDensity = function (t, xn)
local x, y, z = xn[1], xn[2], xn[3]
local sourceFloor = 1e10
if math.abs(z) < Lz/4 then
 near the midplane, the density source is a Gaussian
return math.max(math.exp((xxSource)^2/(2*lambdaSource)^2), sourceFloor)
else
return 1e40
end
end
sourceTemperature = function (t, xn)
local x, y, z = xn[1], xn[2], xn[3]
if math.abs(xxSource) < 3*lambdaSource then
return 80*eV
else
return 30*eV
end
end
This concludes the Preamble. We now have everything we need to initialize
the Gyrokinetic
App. In this input file, the App initialization consists
of 4 sections:

 App initialization

plasmaApp = Plasma.App {

 Common

...

 Species

...

 Fields

...

 Geometry

...
}
 The Common section includes a declaration of parameters that control the
(configuration space) discretization, and time advancement. This first block of
code in Plasma.App
may specify the periodic directions, the MPI
decomposition, and the frequency with which to output certain diagnostics.

 Common

logToFile = true,  will write simulation output log to gksheath_0.log
tEnd = 10e6,  simulation end time [s]
nFrame = 10,  number of output frames for diagnostics
lower = {R  Lx/2, Ly/2, Lz/2},  configuration space domain lower bounds, {x_min, y_min, z_min}
upper = {R + Lx/2, Ly/2, Lz/2},  configuration space domain upper bounds, {x_max, y_max, z_max}
cells = {4, 1, 8},  number of configuration space cells, {nx, ny, nz}
basis = "serendipity",  basis type (only "serendipity" is supported for gyrokinetics)
polyOrder = 1,  polynomial order of basis set (polyOrder = 1 fully supported for gyrokinetics, polyOrder = 2 marginally supported)
timeStepper = "rk3",  timestepping algorithm
cflFrac = 0.4,  fractional modifier for timestep calculation via CFL condition
restartFrameEvery = .2,  restart files will be written after every 20% of simulation
 Specification of periodic directions
 (1based indexing, so xperiodic = 1, yperiodic = 2, etc)
periodicDirs = {2},  Periodic in y only (y = 2nd dimension)
 The Species section sets up the species to be considered in the simulation.
Each species gets its own Lua table, in which one provides the velocityspace domain and discretization of the species, initial conditions, sources, collisions, boundary conditions, and diagnostics.
In this input file, we initialize gyrokinetic electron and ion species. Since this section is the most involved part of the input file, we will discuss various parts in detail below.

 Species

 Gyrokinetic electrons
electron = Plasma.Species {
evolve = true,  evolve species?
charge = qe,  species charge
mass = me,  species mass
 Speciesspecific velocity domain
lower = {4*vte, 0},  velocity space domain lower bounds, {vpar_min, mu_min}
upper = {4*vte, 12*me*vte^2/(2*B0)},  velocity space domain upper bounds, {vpar_max, mu_max}
cells = {8, 4},  number of velocity space cells, {nvpar, nmu}
 Initial conditions
init = Plasma.MaxwellianProjection {  initialize a Maxwellian with the specified density and temperature profiles
 density profile
density = function (t, xn)
 The particular functional form of the initial density profile
 comes from a 1D singlefluid analysis (see Shi thesis), which derives
 quasisteadystate initial profiles from the source parameters.
local x, y, z, vpar, mu = xn[1], xn[2], xn[3], xn[4], xn[5]
local Ls = Lz/4
local floor = 0.1
local effectiveSource = math.max(sourceDensity(t,{x,y,0}), floor)
local c_ss = math.sqrt(5/3*sourceTemperature(t,{x,y,0})/mi)
local nPeak = 4*math.sqrt(5)/3/c_ss*Ls*effectiveSource/2
local perturb = 0
if math.abs(z) <= Ls then
return nPeak*(1+math.sqrt(1(z/Ls)^2))/2*(1+perturb)
else
return nPeak/2*(1+perturb)
end
end,
 temperature profile
temperature = function (t, xn)
local x = xn[1]
if math.abs(xxSource) < 3*lambdaSource then
return 50*eV
else
return 20*eV
end
end,
scaleWithSourcePower = true,  when source is scaled to achieve desired power, scale initial density by same factor
},
 Collisions parameters
coll = Plasma.LBOCollisions {  LenardBernstein model collision operator
collideWith = {'electron'},  only include selfcollisions with electrons
frequencies = {nuElc},  use a constant (in space and time) collision freq. (calculated in Preamble)
},
 Source parameters
source = Plasma.Source {  source is a Maxwellian with the specified density and temperature profiles
density = sourceDensity,  use sourceDensity function (defined in Preamble) for density profile
temperature = sourceTemperature,  use sourceTemperature function (defined in Preamble) for temperature profile
power = P_src/2,  sourceDensity will be scaled to achieve desired power
diagnostics = {"intKE"},
},
 Nonperiodic boundary condition specification
bcx = {Plasma.ZeroFluxBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}},
Plasma.ZeroFluxBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}}},  use zeroflux boundary condition in x direction
bcz = {Plasma.SheathBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}},
Plasma.SheathBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}}},  use sheathmodel boundary condition in z direction
 Diagnostics
diagnostics = {"M0", "Upar", "Temp", "intM0", "intM1", "intKE", "intEnergy"},
},
 Gyrokinetic ions
ion = Plasma.Species {
evolve = true,  evolve species?
charge = qi,  species charge
mass = mi,  species mass
 Speciesspecific velocity domain
lower = {4*vti, 0},  velocity space domain lower bounds, {vpar_min, mu_min}
upper = {4*vti, 12*mi*vti^2/(2*B0)},  velocity space domain upper bounds, {vpar_max, mu_max}
cells = {8, 4},  number of velocity space cells, {nvpar, nmu}
 Initial conditions
init = Plasma.MaxwellianProjection {  initialize a Maxwellian with the specified density and temperature profiles
 density profile
density = function (t, xn)
 The particular functional form of the initial density profile
 comes from a 1D singlefluid analysis (see Shi thesis), which derives
 quasisteadystate initial profiles from the source parameters.
local x, y, z, vpar, mu = xn[1], xn[2], xn[3], xn[4], xn[5]
local Ls = Lz/4
local floor = 0.1
local effectiveSource = math.max(sourceDensity(t,{x,y,0}), floor)
local c_ss = math.sqrt(5/3*sourceTemperature(t,{x,y,0})/mi)
local nPeak = 4*math.sqrt(5)/3/c_ss*Ls*effectiveSource/2
local perturb = 0
if math.abs(z) <= Ls then
return nPeak*(1+math.sqrt(1(z/Ls)^2))/2*(1+perturb)
else
return nPeak/2*(1+perturb)
end
end,
 temperature profile
temperature = function (t, xn)
local x = xn[1]
if math.abs(xxSource) < 3*lambdaSource then
return 50*eV
else
return 20*eV
end
end,
scaleWithSourcePower = true,  when source is scaled to achieve desired power, scale initial density by same factor
},
 Collisions parameters
coll = Plasma.LBOCollisions {  LenardBernstein model collision operator
collideWith = {'ion'},  only include selfcollisions with ions
frequencies = {nuIon},  use a constant (in space and time) collision freq. (calculated in Preamble)
},
 Source parameters
source = Plasma.Source {  source is a Maxwellian with the specified density and temperature profiles
density = sourceDensity,  use sourceDensity function (defined in Preamble) for density profile
temperature = sourceTemperature,  use sourceTemperature function (defined in Preamble) for temperature profile
power = P_src/2,  sourceDensity will be scaled to achieve desired power
diagnostics = {"intKE"},
},
 Nonperiodic boundary condition specification
bcx = {Plasma.ZeroFluxBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}},
Plasma.ZeroFluxBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}}},  use zeroflux boundary condition in x direction
bcz = {Plasma.SheathBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}},
Plasma.SheathBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}}},  use sheathmodel boundary condition in z direction
 Diagnostics
diagnostics = {"M0", "Upar", "Temp", "intM0", "intM1", "intKE", "intEnergy"},
},
This simulation also requires a source, which models plasma crossing the separatrix. The next part of the Preamble initializes some source parameters, along with some functions that will be used later to set up the source density and temperature profiles.
 Source parameters
P_SOL = 3.4e6  total SOL power, from experimental heating power [W]
P_src = P_SOL*Ly*Lz/(2*math.pi*R*Lpol)  fraction of total SOL power into flux tube domain [W]
xSource = R  source peak radial location [m]
lambdaSource = 0.005  source radial width [m]
 Source density and temperature profiles.
 Note that source density will be scaled to achieve desired source power.
sourceDensity = function (t, xn)
local x, y, z = xn[1], xn[2], xn[3]
local sourceFloor = 1e10
if math.abs(z) < Lz/4 then
 near the midplane, the density source is a Gaussian
return math.max(math.exp((xxSource)^2/(2*lambdaSource)^2), sourceFloor)
else
return 1e40
end
end
sourceTemperature = function (t, xn)
local x, y, z = xn[1], xn[2], xn[3]
if math.abs(xxSource) < 3*lambdaSource then
return 80*eV
else
return 30*eV
end
end
This concludes the Preamble. We now have everything we need to initialize
the Gyrokinetic
App. In this input file, the App initialization consists
of 4 sections:

 App initialization

plasmaApp = Plasma.App {

 Common

...

 Species

...
scaleWithSourcePower = true,  when source is scaled to achieve desired power, scale initial density by same factor
},
 Collisions parameters
coll = Plasma.LBOCollisions {  LenardBernstein model collision operator
collideWith = {'ion'},  only include selfcollisions with ions
frequencies = {nuIon},  use a constant (in space and time) collision freq. (calculated in Preamble)
},
 Source parameters
source = Plasma.MaxwellianProjection {  source is a Maxwellian with the specified density and temperature profiles
isSource = true,  designate as source
density = sourceDensity,  use sourceDensity function (defined in Preamble) for density profile
temperature = sourceTemperature,  use sourceTemperature function (defined in Preamble) for temperature profile
power = P_src/2,  sourceDensity will be scaled to achieve desired power
},
 Nonperiodic boundary condition specification
bcx = {Plasma.ZeroFluxBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}},
Plasma.ZeroFluxBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}}},  use zeroflux boundary condition in x direction
bcz = {Plasma.SheathBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}},
Plasma.SheathBC{diagnostics={"M0", "Upar", "Energy", "intM0", "intM1", "intKE", "intEnergy"}}},  use sheathmodel boundary condition in z direction
 Diagnostics
diagnostics = {"M0", "Upar", "Temp", "intM0", "intM1", "intKE", "intEnergy"},
},
The initial condition for this problem is given by a Maxwellian. This
is specified using init = Plasma.MaxwellianProjection { ... }
,
which is a table with entries for the density and temperature profile
functions (we could also specify the driftSpeed profile) to be used
to initialze the Maxwellian. In this simulation, the initial density
profile takes a particular form that comes from a 1D singlefluid
analysis (see [Shi2019]), which derives quasisteadystate initial
profiles from the source parameters.
By default the sources, specified via
source = Plasma.Source { ... }
, also take the form of Maxwellians.
For the density and temperature profile functions, we use the
sourceDensity
and sourceTemperature
functions defined in the
Preamble. We also specify the desired source power. The source density
is then scaled so that the integrated power in the source matches the
desired power. Therefore, sourceDensity only controls the shape of the
source density profile, not the amplitude. Since the initial conditions
are related to the source, we also scale the initial species density
by the same factor as the source via the scaleWithSourcePower = true
flag in the initial conditions.
Selfspecies collisions are included using a LenardBernstein model
collision operator via the coll = Plasma.LBOCollisions { ... }
table.
For more details about collision models and options, see
Collisions.
Nonperiodic boundary conditions are specified via the bcx
and bcz
tables. For this simulation, we use zeroflux boundary conditions in the
\(x\) (radial) direction, and sheathmodel boundary conditions in the
\(z\) (fieldaligned) direction.
Finally, we specify the diagnostics that should be outputted for each species. These consist of various moments and integrated quantities. For more details about available diagnostics, see Gyrokinetic App: Electromagnetic gyrokinetic model for magnetized plasmas.
 The Fields section specifies parameters and options related to the
field solvers for the gyrokinetic potential(s).

 Fields

 Gyrokinetic field(s)
field = Plasma.Field {
evolve = true,  Evolve fields?
isElectromagnetic = false,  use electromagnetic GK by including magnetic vector potential A_parallel?
 Nonperiodic boundary condition specification for electrostatic potential phi
 Dirichlet in x.
phiBcLeft = { T ="D", V = 0.0},
phiBcRight = { T ="D", V = 0.0},
 Periodic in y. 
 No BC required in z (Poisson solve is only in perpendicular x,y directions)
},
 The Geometry section specifies parameters related to the
background magnetic field and other geometry parameters.

 Geometry

 Magnetic geometry
funcField = Plasma.Geometry {
 Background magnetic field profile
 Simple helical (i.e. cylindrical slab) geometry is assumed
bmag = function (t, xn)
local x = xn[1]
return B0*R/x
end,
quantity is controlled by the nFrame
parameter in the input file.
We can use the Gkeyll postprocessing tool (postgkyl) to visualize the outputs.
Electron density and distribution function¶
First, let’s examine the initial conditions, which are given in output file
sending in _0.bp
. The initial electron density \(n_e(x,y,z)\) is
found in gksheath_electron_M0_0.bp
, where M0
is the label for
the density moment. Let’s look at this file as a function of the \(x\)
and \(z\) coordintes by taking a lineout at \(y=0\) via
pgkyl gksheath_electron_M0_0.bp interp sel z1 0. pl x '$x$' y '$z$'
where we have used the interp
(interpolate)
command to interpolate the DG data onto the grid, and the sel z1 0.
(select) command to make the lineout at \(y=0\)
(z1
refers to the \(y\) coordinate here). The resulting plot looks like
We ran this simulation for 10 \(\mu\text{s}\), and since nframe=10
we have an output frame for each \(\mu\text{s}\) of the simulation.
Let’s look at the final state now, at \(t=10\mu\text{s}\).
pgkyl gksheath_electron_M0_10.bp interp sel z1 0. pl x '$x$' y '$z$'
gives
Seeing as we have run a kinetic calculation, we may wish to examine the velocityspace structure of the distribution function. From postgkyl’s point of view distribution functions are just another dataset, albeit a higher dimensional one. Since we can only produce 1D and 2D plots at the moment we have to select at least 3 of the 5 coordinates at specific values. We will make a 2D plot of velocityspace at t=0 by selecting (x,y,z)=(1.,0,0), which is near the center of the domain, with the following command:
pgkyl gksheath_electron_0.bp interp sel z0 1. z1 0. z2 0. pl x '$v_\parallel$' y '$\mu$' clabel '$f_e(x=1,y=0,z=0,v_\parallel,\mu,t=0)$'
This plot shows that the initial \(f_e\) is Maxwellian. In this example the distribution remains essentially Maxwellian throughout time, so if we were to plot the last frame we would obtain a similar picture.
Sheath potential¶
Now let’s look at the electrostatic potential, \(\phi\). We’d like to see if the sheath potential formed selfconsistently due to our conductingsheath boundary conditions. Let’s look at \(\phi\) along the field line (i.e. along the \(z\) coordinate) by taking lineouts at \(x=1.0\) and \(y=0\).
pgkyl gksheath_phi_10.bp interp sel z0 1. z1 0. pl x '$z$'
gives
Indeed, at the domain ends in \(z\), we have a sheath potential \(\phi_{sh} = 90 \text{ V}\).
We can also make an animation of the evolution of the sheath potential via
pgkyl "gksheath_phi_[09]*.bp" interp sel z0 1. z1 0. anim x '$z$'
Divertor Fluxes¶
pgkyl gksheath_ion_bcZlower_flux_M0_10.bp interp ev 'f[0] 1,2 avg' pl x '$x$'
Here we use ev
to average in the \(y\) and \(z\) direction
(for boundary fluxes, an average in the boundary direction is always
required). This results in
The ion energy (heat) flux profile can similarly be plotted via
pgkyl gksheath_ion_bcZlower_flux_Energy_10.bp interp ev 'f[0] 1,2 avg' pl x '$x$'
Suppose instead of the instantaneous flux, we want the timeaveraged flux over some period of time, perhaps from 510 \(\mu\text{s}\). To compute this, we can use
pgkyl "gksheath_ion_bcZlower_flux_Energy_*.bp" interp collect \
sel z0 5:10 ev 'f[0] 0,2,3 avg' pl x '$x$'
This uses the collect command to aggregate the
frames into a time dimension, which becomes coordinate 0. We then use
sel z0 5:10
to select frames 510. Then we use ev 'f[1] 0,2,3 avg'
to average the data in the 0th (time), 2nd (\(y\)), and 3rd (\(z\))
dimensions. This gives
References¶
[Shi2019]  Shi, E. L., Hammett, G. W., StoltzfusDueck, T., & Hakim, A. (2019). “Fullf gyrokinetic simulation of turbulence in a helical openfieldline plasma”, Physics of Plasmas, 26, 012307. https://doi.org/10.1063/1.5074179 
Fluid example¶
Contents
Physics model and initial conditions¶
This simulation models 2d fast magnetic reconnection in a twofluid (electronion) plasma (see [Birn2001] for the classical paper on this topic and [Wang2015] for our earlier application). The electron and ion fluid spcies are evolved using the fivemoment model
where \(\rho_{s}\), \(\mathbf{u}_{s}\), \(p_{s}\), and \(\mathcal{E}=\frac{3}{2}p_{s}+\frac{1}{2}\rho_{s}\mathbf{u}_{s}^{2}\) are the mass density, velocity, thermal pressure, and total energy of the species \(s\). The electromagnetic fields are evolved and coupled to the plasma species by the usual Maxwell’s equations
The algorithms are summarized in more detail in Moments App: MultifluidmomentMaxwell model, as well as in [Hakim2006], [Hakim2008], and [Wang2020].
The simulation input file is adapted from those used in our earlier work presented in [Wang2015], though only a fivemoment version is presented here. We employ a 2D simulation domain that is periodic in \(x\) (\(L_{x}/2<x<L_{x}/2\)) and is bounded by conducting walls in \(y\) at \(y=\pm L_{y}/2\). Here \(L_{x}=100d_{i0}\), \(L_{y}=50d_{i0}\), where \(d_{i0}=\sqrt{m_{i}/\mu_{0}n_{0}\lefte\right^{2}}\) is the ion inertia length due to \(n_{0}\). The initial equilibrium is a single Harris sheet where magnetic field and number densities are specified by
and
respectively. The total current, \(\mathbf{J}_{0}=\nabla\times\mathbf{B}_{0}/\mu_{0}\), is decomposed according to \(J_{ze0}/J_{zi0}=T_{i0}/T_{e0}\), where the initial temperatures \(T_{e0}\) and \(T_{i0}\) are constant. A sinusoid perturbation is then applied on the magnetic field according to \(\delta\mathbf{B}=\hat{\mathbf{z}}\times\nabla\psi\), where
Input file¶
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190  
 App dependencies

 Load the App and equations to be used
local Moments = require("App.PlasmaOnCartGrid").Moments()
local Euler = require "Eq.Euler"

 Preamble

 Basic constants
local gasGamma = 5./3.  gas gamma
local elcCharge = 1.0  electron charge
local ionCharge = 1.0  ion charge
local elcMass = 1.0  electron mass
local lightSpeed = 1.0  light speed
local mu0 = 1.0  mu_0
 Problemspecific characteristic parameters
local n0 = 1.0  characteristic density variation
local nb_n0 = 0.2  background density vs. varation density
local plasmaBeta0 = 1.0  thermal pressure vs. magnetic pressure
local lambda_di0 = 0.5  transition layer thickness over ion inertial length
local Ti0_Te0 = 5.0  ion temperature vs. electron temperature
local massRatio = 25.0  ion mass vs. electron mass
local vAe0_lightSpeed = 0.5  electron Alfven speed vs. c
local perturbationLevel = 0.1  relative perturbation magnitude
 Derived parameters for problem setup or for diagnostic information
local epsilon0 = 1.0 / lightSpeed^2 / mu0  epsilon_0
local ionMass = elcMass * massRatio  ion mass
local vAe0 = lightSpeed * vAe0_lightSpeed  electron Alfven speed
local vAlf0 = vAe0 * math.sqrt(elcMass / ionMass)  plasma Alfven speed
local B0 = vAlf0 * math.sqrt(n0 * ionMass)  background B field
local OmegaCi0 = ionCharge * B0 / ionMass  ion cyclotron frequency
local OmegaPe0 = math.sqrt(n0 * elcCharge^2 / (epsilon0 * elcMass))
 electron plasma frequency
local de0 = lightSpeed / OmegaPe0  electron inertial length
local OmegaPi0 = math.sqrt(n0 * ionCharge^2 / (epsilon0 * ionMass))
 ion plasma frequency
local di0 = lightSpeed / OmegaPi0  ion inertial length
local lambda = lambda_di0 * di0  transition layer thickness
local psi0 = perturbationLevel * B0  potential of perturbation B field
 Domain and time control
local Lx, Ly = 25.6 * di0, 12.8 * di0  domain lengths
local Nx, Ny = 128, 64  grid size
local tEnd = 25.0 / OmegaCi0  end of simulation
local nFrame = 5  number of output frames at t=tEnd

 App construction

local momentApp = Moments.App {

 Common

logToFile = true,  will log to rt5mgem_0.log
tEnd = tEnd,  stop the simulation at t=tEnd
nFrame = nFrame,  number of output frames at t=tEnd
lower = {Lx/2, Ly/2},  lower limit of domain in each direction
upper = {Lx/2, Ly/2},  upper limit of domain in each direction
cells = {Nx, Ny},  number of cells in each direction
timeStepper = "fvDimSplit",  always use fvDimSplit for fluid simulations,
 i.e., finitevolume dimensionalsplitting
periodicDirs = {1},  periodic boundary condition directions

 Species

 electrons
elc = Moments.Species {
charge = elcCharge, mass = elcMass,  charge and mass
equation = Euler { gasGamma = gasGamma },
 default equation to be evolved
 using 2ndorder scheme
equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" },
 backup equation to be solved to
 retake the time step in case of
 negative density/pressure
init = function (t, xn)  initial condition to set the
 moments of this species
local x, y = xn[1], xn[2]
local TeFrac = 1.0 / (1.0 + Ti0_Te0)
local sech2 = (1.0 / math.cosh(y / lambda))^2
local n = n0 * (sech2 + nb_n0)
local Jz = (B0 / lambda) * sech2
local Ttotal = plasmaBeta0 * (B0 * B0) / 2.0 /n0
local rho_e = n * elcMass
local rhovx_e = 0.0
local rhovy_e = 0.0
local rhovz_e = (elcMass / elcCharge) * Jz * TeFrac
local thermalEnergy_e = n * Ttotal * TeFrac / (gasGamma1.0)
local kineticEnergy_e = 0.5 * rhovz_e * rhovz_e / rho_e
local e_e = thermalEnergy_e + kineticEnergy_e
return rho_e, rhovx_e, rhovy_e, rhovz_e, e_e
end,
evolve = true,  whether to evolve this species
 in the hyperbolic step; it could
 still be evolved in the source
 update step, though
bcy = { Euler.bcWall, Euler.bcWall },  boundary conditions along
 nonperiodic directions; in this
 case, the y direction
},
 ions
ion = Moments.Species {
charge = ionCharge, mass = ionMass,
equation = Euler { gasGamma = gasGamma },
equationInv = Euler { gasGamma = gasGamma, numericalFlux = "lax" },
init = function (t, xn)
local x, y = xn[1], xn[2]
local TiFrac = Ti0_Te0 / (1.0 + Ti0_Te0)
local sech2 = (1.0 / math.cosh(y / lambda))^2
local n = n0 * (sech2 + nb_n0)
local Jz = (B0 / lambda) * sech2
local Ttotal = plasmaBeta0 * (B0 * B0)/ 2.0 / n0
local rho_i = n*ionMass
local rhovx_i = 0.0
local rhovy_i = 0.0
local rhovz_i = (ionMass / ionCharge) * Jz * TiFrac
local thermalEnergy_i = n * Ttotal * TiFrac / (gasGamma  1)
local kineticEnergy_i = 0.5 * rhovz_i * rhovz_i / rho_i
local e_i = thermalEnergy_i + kineticEnergy_i
return rho_i, rhovx_i, rhovy_i, rhovz_i, e_i
end,
evolve = true,
bcy = { Euler.bcWall, Euler.bcWall },
},

 Electromagnetic field

field = Moments.Field {
epsilon0 = epsilon0, mu0 = mu0,
init = function (t, xn)
local x, y = xn[1], xn[2]
local pi = math.pi
local sin = math.sin
local cos = math.cos
local Ex, Ey, Ez = 0.0, 0.0, 0.0
local Bx = B0 * math.tanh(y / lambda)
local By = 0.0
local Bz = 0.0
local dBx = psi0 *(pi / Ly) * cos(2 * pi * x / Lx) * sin(pi * y / Ly)
local dBy = psi0 * (2 * pi / Lx) * sin(2 * pi * x / Lx) * cos(pi * y / Ly)
local dBz = 0.0
Bx = Bx + dBx
By = By + dBy
Bz = Bz + dBz
return Ex, Ey, Ez, Bx, By, Bz
end,
evolve = true,  whether to evolve the field in the hyperbolic step; even
 if this is set to false, the E field could still be
 evolved in the source update step
bcy = { Moments.Field.bcReflect, Moments.Field.bcReflect },
 boundary conditions along the nonperiodic directions,
 in this case, the y direction.
},

 Source equations that couple the plasma fluids and EM fields

emSource = Moments.CollisionlessEmSource {
species = {"elc", "ion"},  plasma species involved
timeStepper = "direct",  timestepper; one of "timecentered",
 "direct", and "exact"
},
}

 Run the App

momentApp:run()

Running the simulation¶
The input file rt5mgem.lua
can be run using the gkyl executable
mpirun n 4 ~/gkylsoft/gkyl/bin/gkyl inputFiles/rt5mgem.lua
assuming gkyl
has been installed in the user’s home directory.
When running this simulation, a user should see the following output
Fri Sep 18 2020 15:53:11.000000000
Gkyl built with e0fd133198bd+
Gkyl built on Sep 11 2020 17:59:49
Initializing PlasmaOnCartGrid simulation ...
Using CFL number 1
Initializing completed in 0.0723601 sec
Starting main loop of PlasmaOnCartGrid simulation ...
** Time step too large! Aborting this step! ** Time step 1250 too large! Will retake with dt 1
Step 0 at time 0. Time step 1. Completed 0%
0123456789 Step 125 at time 125. Time step 1. Completed 10%
0123456789 Step 250 at time 250. Time step 1. Completed 20%
0123456789 Step 375 at time 375. Time step 1. Completed 30%
0123
Postprocessing¶
The output of this simulation is the following set of files:
 Electron fluid moments:
rt5mgem_elc_#.bp
.  Ion fluid moments:
rt5mgem_elc_#.bp
.  Electromagnetic fields:
rt5mgem_field_#.bp
.  Field energy:
rt5mgem_fieldEnergy.bp
.  Diagnostic integrated moments:
rt5mgem_elc_intMom.bp
andrt5mgem_in_intMom.bp
.
Snapshots (frames) are labeled by the number #
at the end of the file name, while integrated diagnostics that are computed as a timeseries, such as the field energy, are written out as a single file.
Since nFrame=1
in the input file, the only frames that are output are 0
, corresponding to the initial condition, and 1
, corresponding to the end of the simulation.
Plotting all physical quantities in a frame¶
The following command plot all five moment term in the 5th electron output frame. Note that the five moment terms are ordered as \(\rho,\,\rho v_x\, \rho v_y, \rho v_z, \mathcal{E}\).
pgkyl rt5mgem_elc_5.bp plot fixaspect
Plotting one physical quantity¶
The following command plot the 3rd moment term in the 5th electron output frame., The 3rd moment is the outofplane momentum, which differs from the electron outofplane current by an coefficient \(q_e/m_e\) .
pgkyl rt5mgem_elc_5.bp c 3 plot fixaspect
Plotting integrated quantities over time¶
Finally, since we performed this simulation at two different resolutions, and interesting diagnostic to look at is a comparison of integrated quantities between the two simulations.
For ease of plotting we have moved the data from the two simulations to two different folders, res1
(\(16^2 \times 16^2\)) and res2
(\(32^2 \times 32^2\)).
Here, we are being agnostic on what a user might have named these two different simulations and labeling them ourselves with postgkyl.
pgkyl rt5mgem_fieldEnergy.bp select c 5 plot logy
References¶
[Hakim2006]  Hakim, A., Loverich, J., & Shumlak, U. (2006). A high resolution wave propagation scheme for ideal TwoFluid plasma equations. Journal of Computational Physics, 219, 418–442. https://doi.org/10.1016/j.jcp.2006.03.036 
[Hakim2008]  Hakim, A. H. (2008). Extended MHD modeling with the tenmoment equations. Journal of Fusion Energy, 27, 36–43. https://doi.org/10.1007/s108940079116z 
[Wang2020]  Wang, L., Hakim, A. H., Ng, J., Dong, C., & Germaschewski, K. (2020). Exact and locally implicit source term solvers for multifluidMaxwell systems. Journal of Computational Physics. https://doi.org/10.1016/j.jcp.2020.109510 
[Birn2001]  Birn, L., et al. (2001). Geospace Environmental Modeling ({GEM}) magnetic reconnection challenge. Journal of Geophysical Research, 106, 3715. 
[Wang2015]  (1, 2) Wang, L., Hakim, A. H. A. H., Bhattacharjee, A., & Germaschewski, K. (2015). Comparison of multifluid moment models with particleincell simulations of collisionless magnetic reconnection. Physics of Plasmas, 22(1), 012108. https://doi.org/10.1063/1.4906063 
gkyl Reference¶
Below is the (hopefully) complete reference to the gkyl simulation tool, whose data can be postprocessed with postgkyl. For instructions on installing postgkyl see gkyl install.
Using gkyl¶
Contents
We now cover the basics of runing gkyl on desktops, clusters, and GPUs. Gkyl can also be used be used to run any Lua scripts, as well as tools provided within gkyl. We also comment on some useful tools provided by ADIOS.
Additional details on the contents of the input files can be found in the Input file file basics page and the pages for the Vlasov, Gyrokinetic and Moment Apps.
The installation placed the gkyl executable in
<INSTALLDIR>/gkylsoft/gkyl/bin/
(where the default <INSTALLDIR>
is home, ~
),
so one typically needs to call gkyl as <INSTALLDIR>/gkylsoft/gkyl/bin/gkyl
. However
most users create an alias so one
can simply call gkyl
. The documentation assumes such alias unless specified otherwise.
The gkyl
command has a builtin help menu. Access it with
gkyl h
Run simulations¶
There are three ways of running simulations with gkyl:
 Serial: using a single core/processes/CPU.
 Parallel: running a multicore simulation (using MPI).
 GPUs: using graphical processing units (GPUs).
The input file has to have some knowledge of which of these modalities you will use. We provide some examples of each of these below.
Serial simulations¶
Suppose you have the kbm.lua input file for a linear
kinetic ballooning mode (KBM) calculation with gyrokinetics. In the Common section
of the App declaration (i.e. between plasmaApp = Plasma.App {
and
electron = Plasma.Species {
) there are two variables, decompCuts
and useShared
.
The refer to the number of MPI decompositions and the use of MPI shared memory, respectively.
For serial simulations one can remove these from the input file, or useShared
must be set to false
, and decompCuts
must be a table with as many 1’s as
there are configuration space dimensions (three in this case). That’s why the input
file contains:
plasmaApp = Plasma.App {
...
decompCuts = {1, 1, 1},  Cuts in each configuration direction.
useShared = false,  If to use shared memory.
...
}
Then one can run the input file in serial with the simple command:
gkyl kbm.lua
By the time it completes, after 54 seconds on a 2015 MacbookPro, this simulation will produce the following output to screen:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39  Thu Sep 17 2020 22:20:16.000000000
Gkyl built with 1b66bd4a21e5+
Gkyl built on Sep 17 2020 22:20:05
Initializing Gyrokinetic simulation ...
Initializing completed in 12.9906 sec
Starting main loop of Gyrokinetic simulation ...
Step 0 at time 0. Time step 1.11219e08. Completed 0%
0123456789 Step 27 at time 3.00286e07. Time step 1.11215e08. Completed 10%
0123456789 Step 54 at time 6.00559e07. Time step 1.1121e08. Completed 20%
0123456789 Step 80 at time 8.89697e07. Time step 1.11204e08. Completed 30%
0123456789 Step 107 at time 1.18994e06. Time step 1.11197e08. Completed 40%
0123456789 Step 133 at time 1.47904e06. Time step 1.11189e08. Completed 50%
0123456789 Step 160 at time 1.77924e06. Time step 1.11179e08. Completed 60%
0123456789 Step 186 at time 2.06828e06. Time step 1.11165e08. Completed 70%
0123456789 Step 213 at time 2.3684e06. Time step 1.11145e08. Completed 80%
0123456789 Step 239 at time 2.65735e06. Time step 1.11121e08. Completed 90%
0123456789 Step 266 at time 2.94849e06. Time step 2.27109e09. Completed 100%
0
Total number of timesteps 267
Solver took 25.14505 sec (0.094176 s/step) (46.493%)
Solver BCs took 2.14804 sec (0.008045 s/step) ( 3.972%)
Field solver took 0.58969 sec (0.002209 s/step) ( 1.090%)
Field solver BCs took 0.20732 sec (0.000776 s/step) ( 0.383%)
Function field solver took 0.00000 sec (0.000000 s/step) ( 0.000%)
Moment calculations took 18.12544 sec (0.067886 s/step) (33.514%)
Integrated moment calculations took 4.57880 sec (0.017149 s/step) ( 8.466%)
Field energy calculations took 0.03020 sec (0.000113 s/step) ( 0.056%)
Collision solver(s) took 0.00000 sec (0.000000 s/step) ( 0.000%)
Collision moments(s) took 0.00000 sec (0.000000 s/step) ( 0.000%)
Source updaters took 0.00000 sec (0.000000 s/step) ( 0.000%)
Stepper combine/copy took 1.39611 sec (0.005229 s/step) ( 2.581%)
Time spent in barrier function 0.14791 sec (0.000554 s/step) ( 0.273%)
[Unaccounted for] 1.86320 sec (0.006978 s/step) ( 3.445%)
Main loop completed in 54.08386 sec (0.202561 s/step) ( 100%)
Thu Sep 17 2020 22:21:23.000000000

These simulation logs contain the following:
Line 1:  start date and time. 
Lines 23:  gkyl repository revision with which this simulation was run, and the date on which the executable was built. 
Line 9:  report the initial time step number, time and initial time step size. 
Lines 1019:  report progress every 1% of the simulation (first column). Then, every 10% of the simulation time, give the number of time steps taken so far, simulation time transcurred, and the latest time step size. 
Lines 2137:  give various metrics regarding the timesteps and wallclock time taken by the simulation, and the time spent on various parts of the calculation. 
Line 39:  Date and time when the simulation finished. 
Also, by default gkyl produces a log file with the format <inputfilename>_0.log
.
If you wish to disable this set logToFile = false,
in the Common section of the App.
Parallel simulation¶
For large problems running on a single CPU can lead to impractical runtimes. In those cases one benefits from parallelizing the simulation over many CPUs. This is accomplished in gkyl by decomposing the (phase) space into MPI domains. Therefore, in order to run parallel simulations you must have a parallel installation of gkyl, as most installations typically are.
Suppose one wishes to run the kinetic ballooning mode (KBM) calculation in
the previous section on a node with 16 cores,
using 4 MPI processes along \(y\) and 4 along \(z\). In this case one must edit the
variable decompCuts
in the Common of the input file to reflect this decomposition:
plasmaApp = Plasma.App {
...
decompCuts = {1, 4, 4},  Cuts in each configuration direction.
useShared = false,  If to use shared memory.
...
}
Once decompCuts
and the rest of the input file is set appropriately, you can run
the simulation with the MPI executable provided by your cluster or MPI implementation
(e.g. mpirun, mpiexec, srun, ibrun). For example, with mpirun we would run the simulation as
mpirun n 16 gkyl kbm.lua
The argument following n
is the total number of MPI processes to launch, in this case
\(4\times4=16\). This clearly requires that your computer/node/job has access to
at least 16 cores.
Note
The number of decompCuts
in any dimension should not exceed the number of cells in that dimension.
Note
 (This feature may be superseeded soon) One can request additional
parallelism in velocity space for kinetic simulations by setting
useShared = true
. This enables MPI shared memory. In this case thedecompCuts
must specify the number of nodes and not number of processors. That is, the total number of processors will be determined fromdecompCuts
and the number of threads per node.
On many computer clusters where one may run parallel simulations one must submit scripts in order to submit a job. This jobscript causes the simulation to be queued so that it runs once resources (i.e. cores, nodes) become available. When resources are finally available the simulation runs in a compute node (instead of the login node).
Jobscripts for some machines are provided below. Note that the installation instructions point to machine scripts for building gkyl on each of these computers. If you need assistance with setting up gkyl in a new cluster, see this or feel free to contact the developers.
Sample submit scripts:
Running on GPUs¶
Gkyl is also capable of running on graphical processing units (GPUs) with minimal modifiation
of an input file that you would use to run on CPUs. Our implementation of GPU capabilities uses
CUDA. At the moment, if gkyl was built with CUDA and the node one is performing the computation
in has a GPU, it will default to running the calculation in a GPU. So given an input file
cudaInputFile.lua
, we would simply run it with
gkyl cudaInputFile.lua
On clusters is often common to submit scripts that queue the job for running on compute nodes (when the resources become available). In fact this is often preferable to sshing into a node if that is even possible. Some sample job scripts for running parallel (CPU) jobs were given in the previous section, and below we provide some sample jobscripts for submitting GPU jobs:
Some usage and development notes regarding gkyl’s GPU capabilities can be found in this repository.
Restarts¶
Sometimes a simulations ends prematurely (e.g. your job’s wallclock time allocation ran out), or perhaps it ended successfully but now you wish to run it longer. In these cases one can restart the simulation.
The first simulation prints out a number of restart files, those ending in _restart.bp
. In
order to begin a second simulation from where the first left off, check the tEnd
and nFrame
variables in the input file. These are defined as absolute times/number of frames, that is, they
specify the final simulation time and number of ouput frames from the beginning of the first
simulation, not relative to the previous simulation.
So suppose we run simulation 1 with the following in the App’s Common section:
momentApp = Moments.App {
...
tEnd = 10.0,
nFrame = 100,
...
}
There are two restart scenarios:
 If the simulation completes successfully, one must increase
tEnd
andnFrame
in order to run the second, restart simulation. Otherwise it will just initialize, realize it does not need to advance any further, and terminate. The first simulation ended prematurely, so
tEnd=10.0
was not reached. One can restart the simulation with the sametEnd
andnFrame
and it will simply try to get there this second time. Or one can increasetEnd
andnFrame
so the second simulation goes farther than the first one intended to.
Once you’ve made the appropriate edits to the input file the second, restart simulation is run by simply appending the word restart after the input file, like
This second, restart simulation will use the _restart.bp
files of the first simulation to
construct an initial condition. Note that it will look for the restart files in the same
directory in which the restart simulation is being run, so typically we run restarts in the same
directory as the first simulation.
Using the fromFile
option¶
The fromFile
option can be used to read data from a file on initialization. This can be used
for initial conditions, sources, and geometry data. The file to be read must have the same prefix
as the input file but can otherwise be named as desired, including the extension (it might be useful
to use a different extension, such as .read
, to avoid accidentally deleting needed files if one
does rm *.bp
).
Handy perks¶
Run Lua with gkyl¶
One can use gkyl to run (almost?) any Lua code. Say for example I find code in the interverse which promises to compute the factors of “Life, the Universe, and Everything” (who wouldn’t want that?). We can take such code, put it in an input file named factors.lua and run it with
gkyl factors.lua
Try it! It’s free!
gkyl Tools¶
A number of additional tools that users and developers may find useful as part of their (Gkeyll) workflow are shipped as gkyl Tools. One such tool, for example, allows us to compare BP (ADIOS) files.
Suppose you ran the plasma beach simulation with the Moment App, using the momBeach.lua input file which contains a variable
local J0 = 1.0e12  Amps/m^3.
in the collisionless electromagnetic source. Let’s assume you were scanning this variable, so
you may choose to create another input file momBeachS.lua which
increases J0
to
local J0 = 1.0e10  Amps/m^3.
If after running momBeachS you are not sure if the results changed at all, you can use the
comparefiles
tool. For example, compare the electromagnetic fields produced at the end of
both simulations with the following command:
gkyl comparefiles a momBeach_field_100.bp b momBeachS_field_100.bp
In this particular example the tool would then print the following to screen:
Checking attr numCells in momBeach_field_100.bp momBeach_field_100s.bp ...
... comparing numCells
Checking attr lowerBounds in momBeach_field_100.bp momBeach_field_100s.bp ...
... comparing lowerBounds
Checking attr upperBounds in momBeach_field_100.bp momBeach_field_100s.bp ...
... comparing upperBounds
Checking attr basisType in momBeach_field_100.bp momBeach_field_100s.bp ...
... comparing basisType
Checking attr polyOrder in momBeach_field_100.bp momBeach_field_100s.bp ...
... comparing polyOrder
Files are different!
So we know that increasing J0
by a factor of a 100 did change the simulation.
Additional documentation of these tools is found in the gkyl Tools reference.
ADIOS tools¶
ADIOS has two handy tools that one may use to explore data files produced by a gkyl
simulation. These are bpls
and bpdump
. We give a brief example of each here, and
expanded descriptions of their capabilities can be found in the
ADIOS documentation
, or using the
bpls h
and bpdump h
commands.
Note that these tools are complimentary to postgkyl’s info command.
bpls¶
bpls
provides a simple view of the structure and contents of a .bp
file. For example,
in the previous section we discussed a 5moment calculation of the
plasma beach problem. Such simulation
produced the file momBeach_field_1.bp
. We can explore this file with
bpls momBeach_field_1.bp
which outputs
double time scalar
integer frame scalar
double CartGridField {400, 8}
It tells us that this file contains three variables, the simulation time
at which this snapshot
was produced, the frame
number, and a Cartesian grid field (CartGridField) for 400 cells which
contains 8 electromagnetic components (3 for electric field, 3 for magnetic field, and the other 2
are used in gkyl’s algorithms). One may dump one of these variables with the additional d
flag.
So if we wish to know the simulation time of this frame, we would use
bpls momBeach_field_1.bp time d
and see it output
double time scalar
5.1e11
Note that for large variables (e.g. CartGridField) dumping can overwhelm the terminal/screen. One
can also slice the dataset and only dump part of it, see bpls h
.
There are also a number of attributes (smaller pieces of timeconstant data), which one can see with
the a
flag:
ws:dir jill$ bpls momBeach_field_1.bp a
double time scalar
integer frame scalar
double CartGridField {400, 8}
string changeset attr
string builddate attr
string type attr
string grid attr
integer numCells attr
double lowerBounds attr
double upperBounds attr
string basisType attr
integer polyOrder attr
string inputfile attr
and you can peek the value of an attribute with bpls <filename> a <attributename> d
.
gkyl Apps¶
Gkeyll Apps are toplevel objects that solve a class of problems. Apps make it easy to setup a problem as the steps in the algorithm are prepackaged. Although the apps are usually flexible enough to support most use cases, there may be problems for which a more finegrained control over the simulation is required. In this case, the user has the option of directly instantiating various Gkeyll objects (grids, fields etc) and calling updaters in a custom handcoded timestepping loop [1]. This process can be complicated and is described elsewhere.
The following apps are available in Gkeyll.
Input file basics¶
Input file structure¶
Although this has been covered in the quickstart, we remind the reader that most input files have a similar structure regardless of which model (App) is being used.
On the whole, input files consist of four parts: the App dependencies, the Preamble, the App initialization, and the App run. So most input files look like

 App dependencies.
...

 Preamble.
...

 App initialization.
...

 App run.
...
We expand on the description of each part in the following sections:
App dependencies
This section loads the App one wishes to use for the simulation. At the moment these can be one of VlasovMaxwell, Gyrokinetic or Moments app:local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell()  Load the Vlasov App.
local Plasma = require("App.PlasmaOnCartGrid").Gyrokinetic()  Load the Gyrokinetic App.
local Plasma = require("App.PlasmaOnCartGrid").Moments()  Load the Moments App.
local Constants = require "Lib.Constants"  Load universal physical Constants.
Preamble
In the Preamble one can declare local variables, functions and any object allowed by Lua which may help the user set up the calculation. Often this section is used to define user input parameters, calculate quantities derived from these inputs and create functions which we may later pass to the App. So for example, the simple Landau damping calculation in the first quickstart contained the following Preamble:permitt = 1.0  Permittivity of free space.
permeab = 1.0  Permeability of free space.
eV = 1.0  Elementary charge, or JouleeV conversion factor.
elcMass = 1.0  Electron mass.
ionMass = 1.0  Ion mass.
nElc = 1.0  Electron number density.
nIon = nElc  Ion number density.
Te = 1.0  Electron temperature.
Ti = Te  Ion temperature.
vtElc = math.sqrt(eV*Te/elcMass)  Electron thermal speed.
vtIon = math.sqrt(eV*Ti/ionMass)  Ion thermal speed.
wpe = math.sqrt((eV^2)*nElc/(permitt*elcMass))  Plasma frequency.
lambdaD = vtElc/wpe  Debye length.
 Amplitude and wavenumber of sinusoidal perturbation.
pertA = 1.0e3
pertK = .750/lambdaD
 Maxwellian in (x,vx)space, given the density (denU), bulk flow
 velocity (flowU), mass and temperature (temp).
local function maxwellian1D(x, vx, den, flowU, mass, temp)
local v2 = (vx  flowU)^2
local vtSq = temp/mass
return (den/math.sqrt(2*math.pi*vtSq))*math.exp(v2/(2*vtSq))
end
App initialization
After loading the App, and setting up the App preliminaries in the Preabmle, we must initialize the App itself. This is accomplished with a table likelocal plasmaApp = Plasma.App {

 Common.
...

 Species.
...

 Fields.
...

 ExternalFields.
...

 Extras.
...
}
 Common has parameters that are common to all Apps and control some aspects of the simulation, most notably the final simulation time and the frames to ouput.
 Species contains a declaration of each plasma species to be considered (e.g. electrons, hydrogen ions, neutrals).
 Fields specifies the electrostatic or electromagnetic fields to be included in the simulation.
 ExternalFields: Some simulations also require the specification of (possibly timedependent) external fields. For gyrokinetics, parameters and functions pertaining to the magnetic geometry are specified here.
 Extras: there are additional features that some simulations may require.
App run
Gkyl input files conclude with a call to the run method of the App, in order to get the simulation running once the input file is called by the gkyl executable. This App run command looks likeplasmaApp:run()
The input file Common¶
As mentioned in the previous section, the App initialization has a section called the Common, which contains parameters common to all the Apps. Here we describe what these possible entries are and their default value (if a default value is not given it means that the user must provide this parameter).
Parameter  Description  Default 

tEnd  Final simulation time.  
lower  Table with configuration space coordinates of lower boundaries.  
upper  Table with configuration space coordinates of upper boundaries.  
cells  Table with number of cells along configuration space each direction.  
nFrame  Number of frames of data to write. Initial conditions are always written. For more finegrained control over species and field output, see below.  
periodicDirs  Periodic directions. Note: X is 1, Y is 2 and Z is 3.  { } 
basis  Basis functions to use. One of "serendipity" , "tensor" or "maximalorder" . 
"serendipity" 
polyOrder  Polynomial order of the basis.  0 
basis  Basis functions to use. One of "serendipity" , "tensor" or "maximalorder" . 
"serendipity" 
decompCuts  For parallel simulations: Table with number of processors to use in each configuration space direction.  { } 
useShared  For parallel simulations: Set to true to use MPI shared memory. 
false 
maximumDt  Largest time step size allowed.  tEndtStart 
suggestedDt  Initial suggested timestep. Adjusted as simulation progresses.  maximumDt 
cflFrac  Fraction (usually 1.0) to multiply CFL determined timestep.  1.0, or 2.0 for timeStepper = "rk3s4" . 
cfl  CFL number to use in determining the time step. This parameter should be avoided and cflFrac used instead.  cflFrac/(2*polyOrder+1) 
timeStepper  One of "rk1" (first order RungeKutta), "rk2" (SSPRK2), "rk3"
(SSPRK3) or “rk3s4” (SSPRK3 with 4 stages) or "fvDimSplit" . 
"rk3" 
restartFrameEvery  Frequency with which to write restart files, given as a decimal fraction. Default is every 5% (=0.05) of the simulation, or as frequently as frames are outputted (whichever is largest).  max(0.05, 1./nFrame) 
calcIntQuantEvery  Frequency with which to compute volume integrated quantities, given as a decimal fraction.  Every time step. 
groupDiagnostics  The default (false ) is to output every diagnostic to a separate
file, but if groupDiagnostic=true diagnostics are written to a single
file per frame (per diagnostic app). For example, diagnostics defined on a
positionspace grid such as velocity moments are written to a _gridDiagnostics_ file. 
false 
ioMethod  Method to use for file output. One of "MPI" or "POSIX" . When "POSIX"
is selected, each node writes to its own file in a subdirectory. 
"MPI" 
logToFile  If set to true, log messages are written to log file.  true 
Note
 In general, you should not specify
cfl
orcflFrac
, unless either doing tests or explicitly controlling the timestep. The app will determine the timestep automatically.  The
"rk3s4"
timestepper allows taking twice the timestep as"rk2"
and"rk3"
at the cost of an additional RK stage. Hence, with this stepper a speedup of 1.5X can be expected.  (This feature may be superseeded soon) One can request additional
parallelism in velocity space for kinetic simulations by setting
useShared = true
. This enables MPI shared memory. In this case thedecompCuts
must specify the number of nodes and not number of processors. That is, the total number of processors will be determined fromdecompCuts
and the number of threads per node.
VlasovMaxwell App: VlasovMaxwell equations on a Cartesian grid¶
The VlasovMaxwell
app solves the Vlasov equation on a Cartesian
grid.
where \(f_s\) the particle distribution function, \(\mathbf{a}_s\) is particle acceleration and \(C[f_s,f_{s'}]\) are collisions. This app uses a version of the discontinuous Galerkin (DG) scheme to discretize the phasespace advection, and StrongStability preserving RungeKutta timesteppers (SSPRK) to discretize the time derivative. We use a matrix and quadrature free version of the algorithm described in [Juno2018] which should be consulted for details of the numerics the properties of the discrete system.
For neutral particles, the acceleration can be set to zero. For electromagnetic (or electrostatic) problems, the acceleration is due to the Lorentz force:
The electromagnetic fields can be either specified, or determined from Maxwell or Poisson equations.
Contents
Overall structure of app¶
By CDIM we mean configuration space dimension and by VDIM we mean velocity space dimension. NDIM = CDIM+VDIM is the phasespace dimension. Note that we must have VDIM>=CDIM.
The overall structure of the app is as follows
local Vlasov = (require "App.PlasmaOnCartGrid").VlasovMaxwell()
vlasovApp = Vlasov.App {
 basic parameters
 description of each species: names are arbitrary
elc = Vlasov.Species {
 species parameters
},
 fields (optional, can be omitted for neutral particles)
field = Vlasov.Field {
 field parameters
},
}
 run application
vlasovApp:run()
Note that if the app’s run
method is not called, the simulation
will not be run, but the simulation will be initialized and the
initial conditions will be written out to file.
Basic parameters¶
The app takes the following basic parameters. Parameters that have default values can be omitted.
Parameter  Description  Default 

logToFile  If set to true, log messages are written to log file  true 
tEnd  End time of simulation  
suggestedDt  Initial suggested timestep. Adjusted as simulation progresses.  tEnd/nFrame 
nFrame  Number of frames of data to write. Initial conditions are always written. For more finegrained control over species and field output, see below.  
lower  CDIM length table with lowerleft configuration space coordinates  
upper  CDIM length table with upperright configuration space coordinates  
cells  CDIM length table with number of configuration space cells  
basis  Basis functions to use. One of “serendipity” or “maximalorder”  
polyOrder  Basis function polynomial order  
cfl  CFL number to use. This parameter should be avoided and cflFrac used instead.  Determined from cflFrac 
cflFrac  Fraction (usually 1.0) to multiply CFL determined timestep.  Determined from timeStepper 
timeStepper  One of “rk2” (SSPRK2), “rk3” (SSPRK3) or “rk3s4” (SSPRK3 with 4 stages). For the last, cflFrac is 2.0  “rk3” 
ioMethod  Method to use for file output. One of “MPI” or “POSIX”. When “POSIX” is selected, each node writes to its own file in a subdirectory. Depending on your system “MPI_LUSTRE” may be available and, if so, should be preferred.  “MPI” 
decompCuts  CDIM length table with number of processors to use in each configuration space direction.  { } 
useShared  Set to true to use shared memory. 
false 
periodicDirs  Periodic directions. Note: X is 1, Y is 2 and Z is 3.  { } 
field  Type of field solver to use. See details below. This is optional and if not specified no force terms will be evolved, i.e. the particles will be assumed to be neutral.  nil 
speciesname  Species objects. There can be more than one of these. See details below. 
Note
 In general, you should not specify
cfl
orcflFrac
, unless either doing tests or explicitly controlling the timestep. The app will determine the timestep automatically.  When
useShared=true
thedecompCuts
must specify the number of nodes and not number of processors. That is, the total number of processors will be determined fromdecompCuts
and the number of threads per node.  The “rk3s4” timestepper allows taking twice the timestep as “rk2” and “rk3” at the cost of an additional RK stage. Hence, with this stepper a speedup of 1.5X can be expected.
Note that the field object must be called “field”. You can also omit the field object completely. In this case, it will be assumed that you are evolving neutral particles and the acceleration will be set to zero (i.e. \(\mathbf{a}_s = 0\) in the Vlasov equation).
Only one field object (if not omitted) is required. At present, the app supports a EM field evolved with Maxwell equations, or a EM field specified as a timedependent function.
Species parameters¶
The Vlasov app works with arbitrary number of species. Each species is
described using the Vlasov.Species
objects. By default every
species in the app is evolved. However, species evolution can be
turned off by setting the evolve
flag to false
. Species can be
given arbitrary names. As the species names are used to label the
output data files, reasonable names should be used.
elc = Vlasov.Species {
 species parameters
},
Parameter  Description  Default 

nDistFuncFrame  These many distribution function outputs will be written during
simulation. If not specified, toplevel nFrame parameter
will be used 
nFrame from toplevel 
nDiagnosticFrame  These many diagnostics outputs (moments etc) will be written
during simulation. If not specified, toplevel nFrame
parameter will be used 
nFrame from toplevel 
charge  Species charge (ignored for neutral particles)  
mass  Species mass (ignored for neutral particles)  
lower  VDIM length table with lowerleft velocity space coordinates  
upper  VDIM length table with upperright velocity space coordinates  
cells  VDIM length table with number of velocity space cells  
decompCuts  VDIM length table with number of processors to use in each velocity space direction.  { } 
init  Function with signature function(t,xn) that initializes the
species distribution function. This function must return a
single value, \(f(x,v,t=0)\) at xn , which is a NDIM
vector. 

bcx  Length two table with BCs in X direction. See details on BCs below.  { } 
bcy  Length two table with BCs in Y direction. Only needed if CDIM>1  { } 
bcz  Length two table with BCs in Z direction. Only needed if CDIM>2  { } 
evolve  If set to false the species distribution function is not
evolved. In this case, only initial conditions for this species
will be written to file. 
true 
diagnostics  List of moments (and integrated moments) to compute for diagnostics. See below for list of moments supported.  { } 
The supported diagnostic moments are, “M0”, “M1i”, “M2ij”, “M2” and “M3i” defined by
In these diagnostics, the index \(i,j\) run over \(1\ldots VDIM\).
Additionally, Gkeyll can calculate the weak moments: bulk velocity “u” and the square of thermal velocity “vtSq”.
The boundary conditions (if not periodic) are specified with the
bcx
etc. tables. Each table must have exactly two entries, one for
BC on the lower edge and one for the upper edge. The supported values
are
Parameter  Description 

Vlasov.Species.bcAbsorb  All outgoing particles leave the domain, and none reenter. 
Vlasov.Species.bcCopy  A zerogradient BC, approximating an open domain 
Vlasov.Species.bcReflect  Particles are specularly reflected (i.e. billiard ball reflection) 
Note that often “reflection” boundary condition is used to specify a symmetry for particles.
For example, for a 1x simulation, to specify that the left boundary is a reflector, while the right an absorber use:
bcx = { Vlasov.Species.bcReflect, Vlasov.Species.bcAbsorb }
Electromagnetic field parameters¶
The EM field object is used as follows
field = Vlasov.Field {
 field parameters
},
Parameter  Description  Default 

nFrame  These many field outputs will be written during simulation. If
not specified, toplevel nFrame parameter will be used 
nFrame from toplevel 
epsilon0  Vacuum permittivity (\(\epsilon_0\))  
mu0  Vacuum permeability (\(\mu_0\))  
mgnErrorSpeedFactor  Factor specifying speed for magnetic field divergence error correction  0.0 
elcErrorSpeedFactor  Factor specifying speed for electric field divergence error correction  0.0 
hasMagneticField  Flag to indicate if there is a magnetic field  true 
init  Function with signature function(t,xn) that initializes the
field. This function must return 6 values arranged as
\(E_x, E_y, E_z, B_x, B_y, B_z\) at \(t=0\) at xn ,
which is a CDIM vector. 

bcx  Length two table with BCs in X direction. See details on BCs below.  { } 
bcy  Length two table with BCs in Y direction. Only needed if CDIM>1  { } 
bcz  Length two table with BCs in Z direction. Only needed if CDIM>2  { } 
evolve  If set to false the field is not evolved. In this case,
only initial conditions will be written to file. 
true 
Note: When doing an electrostatic problem with no magnetic field,
set the hasMagneticField
to false
. This will choose
specialized solvers that are much faster and can lead to significant
gain in efficiency.
The boundary conditions (if not periodic) are specified with the
bcx
etc. tables. Each table must have exactly two entries, one for
BC on the lower edge and one for the upper edge. The supported values
are
Parameter  Description 

Vlasov.Field.bcCopy  A zerogradient BC, approximating an open domain 
Vlasov.Field.bcReflect  Perfect electrical conductor wall 
Functional field parameters¶
To peform “testparticle” simulation one can specify a timedependent electromagnetic field which does not react to particle currents.
externalField = Vlasov.ExternalField {
 field parameters
},
Parameter  Description  Default 

nFrame  These many field outputs will be written during simulation. If
not specified, toplevel nFrame parameter will be used 
nFrame from toplevel 
emFunc  Function with signature function(t, xn) that specifies
timedependent EM field. It should return six values, in order,
\(E_x, E_y, E_z, B_x, B_y, B_z\). 

evolve  If set to false the field is not evolved. In this case,
only initial conditions will be written to file. 
true 
App output¶
The app will write distribution function for each species and the EM fields at specified time intervals. Depending on input parameters specified to the species and field block, different number of distribution functions, fields and diagnostics (moments, integrated quantities) will be written.
The output format is ADIOS BP files. Say your input file is called “vlasov.lua” and your species are called “elc” and “ion”. Then, the app will write out the following files:
vlasov_elc_N.bp
vlasov_ion_N.bp
vlasov_field_N.bp
Where N
is the frame number (frame 0 is the initial
conditions). Note that if a species or the field is not evolved, then
only initial conditions will be written.
In addition to the above, optionally diagnostic data may also be written. For example, the moments files are named:
vlasov_elc_M0_N.bp
vlasov_ion_M0_N.bp
vlasov_elc_M1i_N.bp
vlasov_ion_M1i_N.bp
etc, depending on the entries in the diagnostics
table for
each species. In addition, integrated moments for each species are
written:
vlasov_elc_intMom_N.bp
This file has the timedependent “M0”, three contributions of kinetic energy and the “M2” (integrated over configuration space) stored in them.
For the field, the electromagnetic energy components \(E_x^2\), \(E_y^2\), \(E_z^2\), \(B_x^2\), \(B_y^2\), and \(B_z^2\) (integrated over configuration space) are stored in the file:
vlasov_fieldEnergy.bp
These can be plotted using postgkyl in the usual way.
Examples¶
 Advection in a potential well (Field not evolved)
 Landau damping of Langmuir waves
 Twostream instability
 Three species electrostatic shock problem. See [Pusztai2018] for full problem description.
 Advection of particles in a constant magnetic field. (Field not evolved)
 Weibel instability in 1x2v. See [Cagas2017] for full problem description.
References¶
[Juno2018]  Juno, J., Hakim, A., TenBarge, J., Shi, E., & Dorland, W.. “Discontinuous Galerkin algorithms for fully kinetic plasmas”, Journal of Computational Physics, 353, 110–147, 2018. http://doi.org/10.1016/j.jcp.2017.10.009 
[Pusztai2018]  I Pusztai, J M TenBarge, A N Csapó, J Juno, A Hakim, L Yi and T Fülöp “Low Machnumber collisionless electrostatic shocks and associated ion acceleration”. Plasma Phys. Control. Fusion 60 035004, 2018. https://doi.org/10.1088/13616587/aaa2cc 
[Cagas2017]  P. Cagas, A. Hakim, W. Scales, and B. Srinivasan, “Nonlinear saturation of the Weibel instability. Physics of Plasmas”, 24 (11), 112116–8, 2017 http://doi.org/10.1063/1.4994682 
Gyrokinetic App: Electromagnetic gyrokinetic model for magnetized plasmas¶
The Gyrokinetic
App solves the (electromagnetic) gyrokinetic system on a
Cartesian grid.
Contents
Overall structure of app¶
To set up a gyrokinetic simulation, we first need to load the Gyrokinetic
App package.
This should be done at the top of the input file, via
local Plasma = (require "App.PlasmaOnCartGrid").Gyrokinetic()
This creates a table Plasma
that loads the gyrokinetic species, fields, etc. packages.
The general structure of the input file is then

 App dependencies.
local Plasma = (require "App.PlasmaOnCartGrid").Gyrokinetic()
...

 Preamble.
...

 App initialization.
plasmaApp = Plasma.App {

 Common
...

 Species
electron = Plasma.Species {
 GkSpecies parameters
...
},
 other species, e.g. ions

 Fields
field = Plasma.Field {
 GkField parameters
...
},

 ExternalFields
extField = Plasma.Geometry {
 GkGeometry parameters
...
},
}

 App run.
plasmaApp:run()
Kinetic simulations may take additional parameters in the input file Common, such as
Parameter  Description  Default 

nDistFuncFrame  These many distribution function outputs will be written during
simulation. If not specified, toplevel nFrame parameter
will be used 
nFrame from toplevel 
Species parameters¶
The Gyrokinetic App works with an arbitrary number of species.
Each species should be declared as

 Species
species_name = Plasma.Species {
 GkSpecies parameters
...
},
The species name (species_name
here) is arbitrary, but will be used for naming in diagnostic files, so names like ion
or electron
are common.
Here we describe all possible parameters used to specify a gyrokinetic species. Parameters that have default values can be omitted. Units are arbitrary, but often SI units are used. In the following, VDIM refers to the velocity space dimension, and CDIM refers to the configuration space dimension. The gyrokinetic app works for 1X1V (CDIM=1, VDIM=1), 1X2V, 2X2V, and 3X2V. The velocity coordinates are \((v_\parallel, \mu)\). See [Shi2017] for details.
Parameter  Description  Default 

charge  Species charge  1.0 
mass  Species mass  1.0 
lower  VDIMlength table with lowerleft velocity space coordinates  
upper  VDIMlength table with upperright velocity space coordinates  
cells  VDIMlength table with number of velocity space cells  
decompCuts  NOT CURRENTLY SUPPORTED, no processor decomposition in velocity space allowed  
init  Specifies how to initialize the species distribution function. Use a Projection plugin (see Projections),
or a function with signature function(t,xn) that return a
single value, \(f(t=0,xn[0],xn[1],...)\), where xn is a NDIM vector. 

evolve  If set to false the species distribution function is not evolved. In this case, only initial conditions for this species will be written to file. 
true 
bcx  Length two table with BCs in X direction. See details on BCs below.  { } 
bcy  Length two table with BCs in Y direction. Only needed if CDIM>1  { } 
bcz  Length two table with BCs in Z direction. Only needed if CDIM>2  { } 
coll  Collisions plugin. See Collisions models in Gkeyll.  
source  Specifies a source that is added to the RHS on every timestep. Use a Projection plugin (see Projections),
or a function with signature function(t,xn) that return a
single value, \(S(t,xn[0],xn[1],...)\), where xn is a NDIM vector. 

diagnostics  List of moments and volume integrated moments to compute. See below for list of moments supported.  { } 
Note
 In general, you should not specify
cfl
orcflFrac
, unless either doing tests or explicitly controlling the timestep. The app will determine the timestep automatically.  When
useShared=true
thedecompCuts
must specify the number of nodes and not number of processors. That is, the total number of processors will be determined fromdecompCuts
and the number of threads per node.  The “rk3s4” timestepper allows taking twice the timestep as “rk2” and “rk3” at the cost of an additional RK stage. Hence, with this stepper a speedup of 1.5X can be expected.
Diagnostics¶
There are speciesspecific diagnostics available, which mainly consist of moments of the distribution function and integrals (over configurationspace) of these moments. There are also additional species diagnostics which serve as metrics of positivity and collisionsrelated errors.
Currently there are four types of diagnostic moments, defined below. Note that in these definitions \(\mathrm{d}\mathbf{w}=\mathrm{d}v_\parallel\) or \(\mathrm{d}\mathbf{w}=(2\pi B_0/m)\mathrm{d}v_\parallel\mathrm{d}\mu\) depending on whether it is a 1V or a 2V simulation. We also use the notation \(d_v\) to signify the number of physical velocityspace dimensions included, i.e. \(d_v=1\) for 1V and \(d_v=3\) for 2V. Also, \(v^2=v_\parallel^2\) for 1V and \(v^2=v_\parallel^2+2\mu B_0/m\) for 2V.
 Velocity moments of the distribution function, written as functions of configurationspace position on each diagnostic frame. The options are
M0
: number density, \(n = M_0 = \int\mathrm{d}\mathbf{w}~f\).M1
: parallel momentum density, \(M_1=\int\mathrm{d}\mathbf{w}~v_\parallel f\).M2
: energy density, \(M_2 = \int\mathrm{d}\mathbf{w}~v^2 f\).Upar
: parallel flow velocity, \(u_\parallel= M_1/n\).Temp
: temperature, \(T = (m/d_v)(M_2  M_1 u_\parallel)/n\)Beta
: plasma beta, \(\beta = 2\mu_0 nT/B^2\)Energy
: particle energy density (kinetic + potential), \(\mathcal{E}_H = \int\mathrm{d}\mathbf{w}~H f\), where \(H = mv^2/2 + q\phi\) is the Hamiltonian.
 Velocity moments integrated over configurationspace, written as timeseries. The options are
intM0
: particle number, \(N = \int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~f\)intM1
: parallel momentum, \(U = \int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~v_\parallel f\)intM2
: \(\int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~v^2 f\)intKE
: kinetic energy, \(\mathcal{E}_K = ({m}/{2})\int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~v^2 f\)intEnergy
: total (kinetic + potential) energy, \(E_H = \int\mathrm{d}\mathbf{x}\mathrm{d}\mathbf{w}~H f\), where \(H = mv^2/2 + q\phi\) is the Hamiltonian.
Boundary flux diagnostics¶
One can request diagnostics of the fluxes through nonperiodic boundaries by providing a
diagnostics = {}
table to the boundary condition Apps. For example, sheath boundary
conditions along z would do this via
bcz = {Plasma.SheathBC{diagnostics={"M0"}}, Plasma.SheathBC{diagnostics={"M0"}}}
in order to request the particle flux through the sheaths. Another example is provided in the gyrokinetic Quickstart page.
The boundary fluxes are computed via integrals of the time rates of change computed in the ghost cells. If we consider a simple phasespace advection equation in 2X2V without any forces
the weak form used by the algorithm is obtained by multiplying this equation by a test function \(\psi\) and integrating over phase space in a single cell. After an integration by parts one obtains
where the hat means that a numerical flux is constructed, and \(\mathrm{d}\mathbf{z}=\mathrm{d}\mathbf{x}\,\mathrm{d}\mathbf{v}\). In ghost cells only the surface terms corresponding to fluxes through the physical domain boundaries are computed. This means tha in the ghost cell at the upper boundary along \(x\), for example
This is phasespace flux through the upper \(x\) boundary during a stage of the PDE solver. For RungeKutta steppers one must form a linear combination of these fluxes from every stage in the same manner as the time rates of change are combined for forward time stepping. For the sake of simplicity here we just assume a single forward Euler step, and define phasespace flux during a single time step through the upper \(x\) boundary as
where the volume factor \(V\) arises from the phasespace integral on the left side of equation (1). Note that these integrals are over a single cell, and that the quantity \(\Gamma_{\mathbf{z},x_+}\) is phasespace field, \(\Gamma_{\mathbf{z},x_+}=\Gamma_{\mathbf{z},x_+}(\mathbf{x},\mathbf{v})\).
With this boundary flux in mind, if one requests the particle density of the boundary flux
through diagnosticBoundaryFluxMoments={GkM0}
the diagnostic would be computed as
This yields the rate of number density crossing the upper \(x\) boundary (per celllength in
the \(x\) direction of the ghost cell). In order to compute
the number of particles per unit time crossing the upper \(x\) boundary
(diagnosticIntegratedBoundaryFluxMoments={intM0}
) we simply integrate the above quantity over
\(y\) (and multiply it by the \(x\)cell length of the ghost cell)
The final detail is that the files created by these diagnostics contain the fluxes through the boundary accumulated since the last snapshot (frame), not since the beginning of the simulation.
References¶
[Shi2017]  Shi, E. L., Hammett, G. W., StolzfusDueck, T., & Hakim, A. (2017). Gyrokinetic continuum simulation of turbulence in a straight openfieldline plasma. Journal of Plasma Physics, 83, 1–27. http://doi.org/10.1017/S002237781700037X 
Contents
Moments App: MultifluidmomentMaxwell model¶
Summary of model equations¶
The Moment
app solves highmoment multifluid equations on a Cartesian grid, coupled to Maxwell’s equations through the Lorentz force.
 Each fluid could be either fivemoment, where the plasma pressure is assumed to a scalar (see [Hakim+2006]), or tenmoment, where the full anisotropic and nongyrotropic plasma pressure tensors (see [Hakim2008]) are evolved.
 The hyperbolic system is solved in converative forms of the fivemoment (Euler) equations, and/or the tenmoment equations, and the perfectly hyperbolic Maxwell’s equations.
 The sources that couples the plasma species momenta and electromagnetic fields are described here and more comprehensively in [Wang+2020].
This App solves the hyperbolic and source parts parts of the coupled system separately and apply high accuracy schemes on both.
Ignoring sources, the homogeneous equations of the fluid moments and electromagneic fields are solved separately using a highresolution wavepropagation finitevolume method described in [Hakim+2006]. The main timestep size constraint comes from the speed of light.
The sources are evolved using a locally implicit, timecentered solver to step over the constraining scales like the Debye length and plasma frequency, etc. See Handling twofluid fivemoment and tenmoment source terms or [Wang+2020] for more details.
Additonial sources can be added if needed, for example, for the tenmoment model, we may apply the closure to relax the pressure tensor towards a scalar pressure (see [Wang+2015]).
We then apply an Strangtype operatorsplitting sequence to combine the hyperbolic and source parts to achive secondorder accuracy in time:
\[\exp\left(\mathcal{L}_{S}\Delta t/2\right)\exp\left(\mathcal{L}_{H}\Delta t\right)\exp\left(\mathcal{L}_{S}\Delta t/2\right).\]Here, we represent the homogeneous update schematically as the operator \(\exp\left(\mathcal{L}_{H}\Delta t\right)\) and the source update as \(\exp\left(\mathcal{L}_{S}\Delta t\right)\).
Overall structure of the Moments app¶

 Preamble 

 The Moments app wraps fluid and field objects, and tells the
 the program how to evolve and couple them
local Moments = require("App.PlasmaOnCartGrid").Moments()
local TenMoment = require "Eq.TenMoment"  TenMoment or Euler
 Create the app
momentApp = Moments.App {

 COMMON 

 basic parameters, e.g., time step, grid, domain decomposition
 Description of each species: names are arbitrary
electron = Moments.Species {
 species parameters, equations, and boundary conditions
},
 Repeat to add more species
hydrogen = Moments.Species { ... },
oxygen = Moments.Species { ... },
 EM fields (optional, can be omitted for neutral fluids)
field = Moments.Field {
 EM field parameters, equations, and boundary conditions
},
 Basic source that couple the fluids and EM fields
emSource = Moments.CollisionlessEmSource {
 names of the species to be coupled
species = {"electron", "hydorgen", "oxygen"},
 other specifications
},
 Additional sources if needed
elc10mSource = Moments.TenMomentRelaxSource {
species = {"elctron"},
 other specifications
},
}
 run the app
momentApp:run()
Examples¶
 Fivemoment modeling of the GEM challenge magnetic reconnection problem.
 Tenmoment modeling of the GEM challenge magnetic reconnection problem. This simulation uses a simplified closure appropriate for this problem.
Basic parameters¶
Parameter  Description  Default 

logToFile 
If set to true, log messages are written to log file  true 
tEnd 
End time of simulation  
suggestedDt 
Initial suggested timestep. Adjusted as simulation progresses.  tEnd/nFrame 
nFrame 
Number of frames of data to write. Initial conditions are always written. For more finegrained control over species and field output, see below.  
lower 
CDIM length table with lowerleft configuration space coordinates  
upper 
CDIM length table with upperright configuration space coordinates  
cells 
CDIM length table with number of configuration space cells  
cfl 
CFL number to use. This parameter should be avoided and ``cflFrac`` used instead.  Determined from cflFrac 
cflFrac 
Fraction (usually 1.0) to multiply CFL determined timestep.  Determined from timeStepper 
maximumDt 
Hard limit of time step size.  tEndtStart 
timeStepper 
The multifluidMaxwell model currently only supports the dimensional
splitting finitevolume method, i.e., "fvDimSplit" . 
"fvDimSplit" 
decompCuts 
CDIM length table with number of processors to use in each configuration space direction.  { } 
useShared 
Set to true to use shared memory. 
false 
periodicDirs 
Periodic directions. Note: X is 1, Y is 2 and Z is 3. E.g., {2} sets
the Y direction to be periodic. 
{ } 
Note
 In general, you should not specify
cfl
orcflFrac
, unless either doing tests or explicitly controlling the timestep. The app will determine the timestep automatically.  When
useShared=true
thedecompCuts
must specify the number of nodes and not number of processors. That is, the total number of processors will be determined fromdecompCuts
and the number of threads per node.
Species parameters¶
Parameter  Description  Default 

charge 
Species charge (ignored for neutral particles)  
mass 
Species mass (ignored for neutral particles)  
equation 
The type of default moment equation for this species, e.g.,
`Euler {gasGamma=5/3}` , `equation = TenMoment {}` . If domain
invariance is violated, i.e., negative density/pressure occurs, the
step is retaken using the `equationInv` method that is supposed to
guarantee positivity but is more diffusive. 

equationInv 
Backup equation that guarantees positivity in case it is violated when
the default `equation` is used. Examples are:
`equationInv = Euler { gasGamma = gasGamma, numericalFlux = 'lax' }` ,
`equationInv = TenMoment { numericalFlux = "lax" }` . 

init 
Function with signature function(t,xn) that initializes the
species moments. This function return n values, where n is the number
of moments for this species. 

bcx 
Length two table with BCs in X direction. See details on BCs below.  { } 
bcy 
Length two table with BCs in Y direction. Only needed if CDIM>1  { } 
bcz 
Length two table with BCs in Z direction. Only needed if CDIM>2  { } 
evolve 
If set to false the moments are not evolved in the hyperbolic part,
but could be modified in the source updater. In this case, by default
only initial conditions for this species will be written to file. To
force writing to file as usual, set the forceWrite option to true. 
true 
forceWrite 
If set to true the moments are written to file even if evolve is
set to false . 
false 
Electromagnetic field parameters¶
Parameter  Description  Default 

nFrame 
These many field outputs will be written during simulation. If
not specified, toplevel nFrame parameter will be used 
nFrame from toplevel 
epsilon0 
Vacuum permittivity (\(\epsilon_0\))  
mu0 
Vacuum permeability (\(\mu_0\))  
mgnErrorSpeedFactor 
Factor specifying speed for magnetic field divergence error correction  0.0 
elcErrorSpeedFactor 
Factor specifying speed for electric field divergence error correction  0.0 
init 
Function with signature function(t,xn) that initializes the
field. This function must return 6 values arranged as
\(E_x, E_y, E_z, B_x, B_y, B_z\) at \(t=0\) at xn ,
which is a CDIM vector. 

bcx 
Length two table with BCs in X direction. See details on BCs below.  { } 
bcy 
Length two table with BCs in Y direction. Only needed if CDIM>1  { } 
bcz 
Length two table with BCs in Z direction. Only needed if CDIM>2  { } 
evolve 
If set to false the field is not evolved. In this case,
only initial conditions will be written to file. 
true 
forceWrite 
If set to true the moments are written to file even if evolve is
set to false . 
false 
App output¶
The app will write snapshots of moments for each species and the EM fields at specified time intervals. Diagnostics like integrated fluid moments and field energy are recorded for each timestep and written in one file for each species/field object.
The output format is ADIOS BP files. Say your input file is called “5m.lua” and your species are called “elc” and “ion”. Then, over specified time invertals the app will write out the following files:
5m_elc_N.bp
5m_ion_N.bp
5m_field_N.bp
Where N
is the frame number (frame 0 is the initial
conditions). Note that if a species or the field is not evolved, then
only initial conditions will be written unless the forceWrite
option
is set to true
.
In addition, integrated moments for each species are written:
vlasov_elc_intMom_N.bp
For the field, the electromagnetic energy components \(E_x^2\), \(E_y^2\), \(E_z^2\), \(B_x^2\), \(B_y^2\), and \(B_z^2\) (integrated over configuration space) are stored in the file:
vlasov_fieldEnergy_N.bp
These can be plotted using postgkyl in the usual way.
References¶
[Hakim+2006]  (1, 2) Hakim, A., Loverich, J., & Shumlak, U. (2006). A high resolution wave propagation scheme for ideal TwoFluid plasma equations. Journal of Computational Physics, 219, 418–442. https://doi.org/10.1016/j.jcp.2006.03.036 
[Hakim2008]  Hakim, A. H. (2008). Extended MHD modeling with the tenmoment equations. Journal of Fusion Energy, 27, 36–43. https://doi.org/10.1007/s108940079116z 
[Wang+2020]  (1, 2) Wang, L., Hakim, A. H., Ng, J., Dong, C., & Germaschewski, K. (2020). Exact and locally implicit source term solvers for multifluidMaxwell systems. Journal of Computational Physics. https://doi.org/10.1016/j.jcp.2020.109510 
[Wang+2015]  Wang, L., Hakim, A. H. A. H., Bhattacharjee, A., & Germaschewski, K. (2015). Comparison of multifluid moment models with particleincell simulations of collisionless magnetic reconnection. Physics of Plasmas, 22(1), 012108. https://doi.org/10.1063/1.4906063 
Footnotes
[1]  The concept of apps was introduced in Gkeyll 2.0. In the 1.0 version of the code the input files contained the complete simulation cycle. Although Gkeyll can still be used in that fashion, it is best to switch to the app model as it simplifies the user input considerably. 
App plugins¶
There are additional features that may be plugged into various Apps.
Collision models in Gkeyll¶
In Gkeyll we currently have two different collision operators for use in kinetic models: the Bhatnagar–Gross–Krook (BGK) and the Dougherty operators. We referred to the latter as the LBO for the legacy of LenardBernstein. Its implementation in Gkeyll is detailed in [Hakim2020] [Francisquez2020].
Contents
BGK collisions¶
The BGK operator [Gross1956] for the effect of collisions on the distribution of species \(f_s\) is
where the sum is over all the species. The distribution functon \(f_{Msr}\) is the Maxwellian
with the primitive moments \(\mathbf{u}_{sr}\) and \(v_{tsr}^2\) properly defined to preserve some properties (such as conservation), and \(d_v\) is the number of velocityspace dimensions. For selfspecies collisions \(\mathbf{u}_{sr}=\mathbf{u}_s\) and \(v_{tsr}^2=v_{ts}^2\). For multispecies collisions we follow an approach similar to [Greene1973] and define the crossspecies primitive moments as
but contrary to Greene’s definition of \(\alpha_E\), we currently use in Gkeyll the following expression
Little guidance is provided by Greene as to how to choose \(\beta\), although it seems clear that \(1<\beta\). In Gkeyll the default value is \(\beta=0\), but the user can specify it in the input file (explained below). We have introduced the additional quantity \(\delta_s\) (which Greene indirectly assumed to equal 1) defined as
The BGK operator can be used with both the VlasovMaxwell solver and the gyrokinetic solver.
Dougherty collisions¶
The Doughery (LBO) model for collisions [Dougherty1964] in Gkeyll is given by
In this case we compute the crossprimitive moments by a process analogous to Greene’s with the BGK operator, yielding the following formulas for the cross flow velocity and thermal speed:
with \(\alpha_E\) defined in the BGK section above. The LBO used by the gyrokinetic solver is
Collisions in Gkeyll input files¶
Users can specify collisions in input files by adding an additional Lua table within each species one wishes to add collisions to. The collision frequency can be constant, varying in space and time, or it can also have a userdefined profile.
Constant collisionality¶
An example of adding
LBO collisions (for BGK collisions simply replace LBOcollisions
with
BGKCollisions
) to a species named ‘elc’ is
elc = Plasma.Species {
charge = q_e, mass = m_e,
 Velocity space grid.
...
 Initial conditions.
...
evolve = true,
 Collisions.
coll = Plasma.LBOCollisions {
collideWith = { "elc" },
frequencies = { nu_ee },
},
},
If there were another species, say one named ‘ion’, this ‘elc’ species could
be made to collide with ‘ion’ by adding ‘ion’ to the collideWidth
table:
coll = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
frequencies = { nu_ee, nu_ei },
},
The constant collision frequencies nu_ee
and nu_ei
need to be previously
computed/specified in the input file. The user can specify the value of \(\beta\)
in the above formulas for the crossspecies primitive moments (\(\mathbf{u}_{sr}\)
and \(v_{tsr}^2\)) by specifying the variable betaGreene
in the collisions
table (if the user does not specify it, betaGreene=0.0
is assumed) like
coll = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
frequencies = { nu_ee, nu_ei },
betaGreene = 0.9
},
In some cases the user may be interested in colliding species ‘elc’ with species ‘ion’, but not collide species ‘ion’ with species ‘elc’. Gkeyll supports this combination, but since the formulas for crossspecies primitive moments involve both \(\nu_{ei}\) and \(\nu_{ie}\), the code will default to assuming \(\nu_{ie}=m_e\nu_{ei}/m_i\). Note however that this scenario is not energy conserving: for exact energy conservation, one must include the effect of binary collisions on both species.
It is also possible to specify both LBO and BGK collisions between different binary pairs in a single input file. For example, if there are three species ‘elc’, ‘ion’ and ‘neut’, the ‘elc’ species could be made collide with both ‘ion’ and ‘neut’ as follows:
cColl = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
frequencies = { nu_ee, nu_ei },
},
nColl = Plasma.BGKCollisions {
collideWith = { "neut" },
frequencies = { nu_en },
},
If no collisionality is specified in the input file, it is assumed that the user desires Gkeyll to build a spatiallyvarying collisionality from scratch using a Spitzerlike formula for \(\nu_{sr}\) (explained below).
Spatially varying collisionality¶
Currently there are three ways to run simulations with a spatially varying collisionality. All of these options lead to a spatially varying, cellwise constant collisionality. We will be adding support for variation of the collisionality within a cell in the future.
Option A¶
The simplest way to run with spatially varying collisionality is to not specify
the table frequencies
. In this case the code computes \(\nu_{sr}\)
according to
where \(\nu_{\mathrm{frac}}\) is a scaling factor, the Coulomb logarithm is defined as
and the \(\alpha\)sum is over all the species. For VlasovMaxwell simulations we do not add the correction due to gyromotion (\(\omega_{c\alpha}=0\) here). The relative velocity here is computed as \(u^2=3v_{tr}^2+3v_{ts}^2\), the reduced mass is \(m_{sr} = m_sm_r/\left(m_s+m_r\right)\), and \(\omega_{p\alpha}\) is the plasma frequency computed with the density and mass of species \(\alpha\). Simpler formulas for the Coulomb logarithm can be easily generated by developers if necessary.
The formulas above assume all the plasma quantities and universal constants are in SI units. The user can provide a different value for these variables by passing them to the collisions table in the input files, as shown here:
coll = Plasma.LBOCollisions {
collideWith = { "elc", "ion" },
epsilon0 = 1.0,  Vacuum permitivity.
elemCharge = 1.0,  Elementary charge value.
hBar = 1.0,  Planck's constant h/2pi.
},
Additionally the user can pass the scaling factor \(\nu_{\mathrm{frac}}\) by
specifying nuFrac
in the collisions table.
Option B¶
Another way to use a spatially varying collisionality is to pass a reference
collisionality normalized to a combination of the density and thermal speed of the
colliding species. This normalized collisionality, is defined as
\(\nu_{srN}=\nu_{sr0}\left(v_{ts0}^2+v_{tr0}^2\right)^{3/2}/n_{r0}\) and one
provides through normNu
in the collisions table as shown below:
elc = Plasma.Species {
...
coll = Plasma.LBOCollisions {
collideWith = { "ion" },
normNu = { nu_ei*((vte^2+vti^2)^(3/2))/n_i0 }
},
},
where nu_ei
, vte
, vti
, n_e0
are computed in the Preamble of the
input file and it is up to the user to ensure that these all have consistent units.
Then, in each time step, the collisions will be applied with the following collisionality
Note that if one is using the normNu
feature for selfspecies collisions, one must
still use these formulas. In this case one would specify electronelectron collisions like
elc = Plasma.Species {
...
coll = Plasma.LBOCollisions {
collideWith = { "elc" },
normNu = { nu_ee*((2*(vte^2))^(3/2))/n_e0 }
},
},
Option C¶
The user may also wish to specify their own collisionality profile, so for this purpose
one can pass functions into the frequencies
table in the collisions table.
For example, suppose that one would like to run a simulation with a collisionality that decays exponentially in x. In this case we could create a exponentially decaying function in the preamble and pass it as the collision frequency as follows:
local Plasma = require("App.PlasmaOnCartGrid").VlasovMaxwell
local Constants = require "Lib.Constants"
eps0 = Constants.EPSILON0
eV = Constants.ELEMENTARY_CHARGE
me = Constants.ELECTRON_MASS
n0 = 7e19  Number density [1/m^3].
Te0 = 100*eV  Electron temperature [J].
 Reference electron collision frequency (at x=0).
logLambdaElc = 24.0  0.5*math.log(n0/1e6) + math.log(Te0/eV)
nu_ee = logLambdaElc*(eV^4)*n0
/(12*math.sqrt(2)*(math.pi^(3/2))*(eps0^2)*math.sqrt(me)*(Te0^(3/2)))
local function nu_eeProfile(t, xn)
local x = xn[1]
return nu_ee*math.exp(x)
end
vlasovApp = Plasma.App {
...
elc = Plasma.Species {
...
 Collisions.
coll = Plasma.LBOCollisions {
collideWith = { "elc" },
frequencies = { nu_eeProfile },
}
}
...
}
 Run application.
vlasovApp:run()
At present all the frequencies
must either be constant numbers or functions. We do not
yet support having a combination of the two in the same collisions table.
Comments on stability¶
The are known issues with the implementation of the collision operators in Gkeyll. One of them, for example, is that we do not have a positivy preseving algorithm for the LBO. Positivity issues are often accompanied by large flows or negative temperatures and/or densities. For this reason we have taken three precautions:
 Calculation of primitive moments \(\mathbf{u}_{sr}\) and \(v_{tsr}^2\) is carried out using cellaverage values if the number density is nonpositive at one of the corners of that cell.
 The collision term is turned off locally if the flow velocity \(\mathbf{u}_{sr}\) is greater than the velocity limits of the domain, or if \(v_{tsr}^2\) is negative.
 The collision frequency \(\nu_{sr}\) is locally set to zero if the cellaverage values of \(n_r\) or \(v_{tsr}^2\) are negative.
We track the number of cells in which precaution 2 is used, and for stable simulations this is typically small (a few percent or less). Further discussion of why these precautions are necessary appears in [Hakim2020].
Examples¶
We offer two full examples of the use of collisions. One in VlasovMaxwell and one in Gyrokinetics.
Example 1: 1x1v collisional relaxation¶
Consider an initial distribution function in 1x1v phase space given by a Maxwellian and a large bump in its tail
Suppose we wish to collisionally relax this initial state, without the influence of collisionless terms. That is, we wish to evolve this distribution function according to equation (1). In this case our input file will use the VlasovMaxwell App (for 1x1v it would be equivalent to use the Gyrokinetic App), and we define the distribution in equation (2) in the Preamble via the function
 Maxwellian with a Maxwellian bump in the tail.
local function bumpMaxwell(x,vx,n,u,vth,bN,bU,bVth,bL,bS)
local vSq = ((vxu)/(math.sqrt(2.0)*vth))^2
local vbSq = ((vxbU)/(math.sqrt(2.0)*bVth))^2
return (n/math.sqrt(2.0*math.pi*vth))*math.exp(vSq)
+(bN/math.sqrt(2.0*math.pi*bVth))*math.exp(vbSq)/((vxbL)^2+bS^2)
end
In this case we chose constants for all densities, flow speed and temperatures. We
also set the charge to 0. Under these conditions the collisionless terms have no effect,
but we can explicitly turn them off with the evolveCollisionless
flag. We will also
request the total integrated bulk flow energy (intM2Flow
) and the total thermal
energy (intM2Thermal
) as diagnostics.
plasmaApp = Plasma.App {
tEnd = 80,  End time.
nFrame = 80,  Number of frames to write.
lower = {0.0},  Configuration space lower coordinate.
upper = {1.0},  Configuration space upper coordinate.
cells = {8},  Configuration space cells.
polyOrder = 2,  Polynomial order.
periodicDirs = {1},  Periodic directions.
 Neutral species with a bump in the tail.
bump = Plasma.Species {
charge = 0.0, mass = 1.0,
 Velocity space grid.
lower = {8.0*vt0}, upper = { 8.0*vt0},
cells = {32},
 Initial conditions.
init = function (t, xn)
local x, v = xn[1], xn[2]
return bumpMaxwell(x,v,n0,u0,vt0,nb,ub,vtb,uL,sb)
end,
evolve = true,  Evolve species?
evolveCollisionless = false,  Evolve collisionless terms?
diagnosticIntegratedMoments = { "intM2Flow", "intM2Thermal" },
 Collisions.
coll = Plasma.LBOCollisions {
collideWith = {'bump'},
frequencies = {nu},
},
},
}
We run this input file with the call
gkyl lboRelax.lua
On a 2015 MacBookPro this ran in 1.5 seconds and produced a screen output like this one.
We can start looking at the data by first, for example, making a movie
of the distribution function as function of time with pgkyl
:
pgkyl "lboRelax_bump_[09]*.bp" interp sel z0 0. anim x '$v$' y '$f(x=0,v,t)$'
(note that postgkyl allows abbreviations,
so interp
= interpolate, sel
= select,
anim
= animate) This command produces the movie given below. We can see that from the
initial, bumpintail state the distribution relaxes to a Maxwellian.
The Maxwellian by the way is the analytic steady state of this operator.
Such relaxation should also take place without breaking momentum or
energy conservation. We can examine the evolution of the total energy
in the system by adding intM2Flow
and intM2Thermal
and plotting
it as a function of time. This is achieved in pgkyl
via:
pgkyl lboRelax_bump_intM2Flow.bp lboRelax_bump_intM2Thermal.bp ev 'f[0] f[1] +' pl x 'time' y 'energy'
As we can see in the figure below, and in particular in the \(10^{14}\) scale of it, the total particle energy is conserved very well. The changes in energy over a collisional period are of the order of machine precision.
Example 2: 1x2v collisional Landau damping¶
We now explore the modification of Landau damping by inclusion of Dougherty collisions. Specifically, we will consider ion acoustic waves with adiabatic electrons. This means that the electron number density simply follows
and our gyrokinetic Poisson equation is simply replaced by the quasineutrality
So there is no need to evolve the electron distribution function. In the Gyrokinetic
App we can specify an adiabatic species using Plasma.AdiabaticSpecies
:
plasmaApp = Plasma.App {
...
adiabaticElectron = Plasma.AdiabaticSpecies {
charge = 1.0, mass = mElc,
temp = Te,
 Initial conditions.. use ion background so that background is exactly neutral.
init = function (t, xn)
return nElc
end,
evolve = false,  Evolve species?
},
...
}
This simulation then only needs to solve the electrostatic gyrokinetic equations for ions
and we do so with an initial condition that contains a sinusoidal perturbation (wavenumber \(k=0.5\)) in the ion density:
If the right side of this equation (4) were zero, this ion acoustic wave would damp at the collisionless rate calculated by Landau (well he did electron Langmuir waves). But collisions will change the picture and we wish to numerically find out how.
This simulation is setup in the ionSound.lua input file. This input file calls for discretizing the ion phase space \([\pi/k,\pi/k]\times[6v_t,6v_t]\times[0,m_i(5v_t^2)/(2B_0)]\) using \(64\times128\times16\) cells and a piecewise linear basis. With a collisionality of \(\nu_=0.005\), the simulation ran on a 2015 MacbookPro in 41 minutes, while a collisionality of \(\nu=0.05\) required 1.35 hours. They were run the command
gkyl ionSound.lua
and produced this screen output.
Note that this is really a linear problem, that is, one can sufficiently model it
with a linearized version of equation (4), using \(f_i=f_{i0}+f_{i1}\),
where the fluctuation \(f_{i1}\) is small compared to the equilibrium
(Maxwellian) \(f_{i0}\). Users may wish to output this fluctuation in time: in
order to to this specify the background with the initBackground
table:
ion = Plasma.Species {
...
 Specify background so that we can plot perturbed distribution and moments.
initBackground = {"maxwellian",
density = function (t, xn)
return nIon
end,
temperature = function (t, xn)
return Ti
end,
},
...
},
This will output the fluctuation to a file with the name format
<simulation>_<species>_f1_#.bp
, where #
stands for the frame number. So for
example, in this ionSound.lua
case it creates files named ionSound_ion_f1_#.bp
.
We can plot this fluctuation along \(v_\parallel\) at \($t=5$\) with
pgkyl "ionSound_ion_f1_10.bp" interp sel z0 0.0 z2 0.0 pl x '$v_\parallel$' y '$f_{i1}(x=0,v_\parallel,\mu=0,t=5)$'
(note that postgkyl allows abbreviations,
so interp
= interpolate, sel
= select,
pl
= plot) which produces the following image
Perhaps most valuable to the physics of this simulation is to see a signature of the
decay of the ion acoustic wave. This simulation produced the integrated squared
electrostatic potential, \(\int\mathrm{d}x\,\phi^2\), which we take as a measure
of the wave energy. It is stored in a file with the name format <simulation>_phiSq.bp
.
If we had run two simulations, ionSound.lua with \(\nu=0.005\)
and ionSoundH.lua with \(\nu=0.05\), we could plot both electrostatic
energies in time with the following pgkyl
command:
pgkyl ionSound_phi2.bp l '$\nu=0.005$' ionSoundH_phi2.bp l '$\nu=0.05$' pl logy f0 x 'time' y 'Integrated $\phi^2$'
Notice that we are giving each file a label to use in the plot with the l
flag. Postgkyl
then produces the following figure
We thus see that the wave energy is decaying as a function of time (the envelope of the curve is going down), and that the rate at which this happens decreases with collisionality. That is, for this case increasing collisionality decreased the damping rate. From this curve we can also read the period of the wave, using the spacing between the dips.
References¶
[Gross1956]  E. P. Gross & M. Krook. Model for collision precesses in gases: smallamplitude oscillations of charged twocomponent systems. Physical Review, 102(3), 593–604 (1956). 
[Greene1973]  J. M. Greene. Improved BhatnagarGrossKrook model of electronion collisions. Physics of Fluids, 16(11), 2022–2023 (1973). 
[Dougherty1964]  J. P. Dougherty. Model FokkerPlanck Equation for a Plasma and Its Solution. Physics of Fluids, 7(11), 1788–1799 (1964). 
[Hakim2020]  (1, 2) A. Hakim, et al. (2020). Conservative Discontinuous Galerkin Schemes for Nonlinear FokkerPlanck Collision Operators. Journal of Plasma Physics Vol 86 No. 4, 905860403 (2020), arXiv:1903.08062. 
[Francisquez2020]  M. Francisquez, et al. (2020). Conservative discontinuous Galerkin scheme of a gyroaveraged Dougherty collision operator. Nuclear Fusion 60 No. 9, 096021 (2020), arxiv:2009.06660. 
Neutral models in Gkeyll¶
Simplified models of plasmaneutral interactions available in Gkeyll. The neutral species is always modeled as a kinetic Boltzmann species using the VlasovMaxwell solver, while the plasma species to which it is coupled can be either VlasovMaxwell or gyrokinetic.
Contents
Electronimpact Ionization¶
This process is given by \(e^{} + n \rightarrow i^{+} + 2e^{}  E_{iz}\), where \(E_{iz}\) is the ionization energy and modeled by collision terms on the RHS of the dynamical equations, as in [Wersal2015]:
where \(f_{M,iz} = f_{M,iz}(n_e, \mathbf{u}_n, v^2_{th,iz})\) is a Maxwellian distribution function that accounts for the lowerenergy elelectrons resulting from this process. (See Neutral species and gyrokinetic plasma species coupling for definitions of Maxwellian distribution function on the VlasovMaxwell and gyrokinetic grids.) \(n_e\) is the electron density, \(\mathbf{u}_n\) is the neutral fluid velocity, and \(v^2_{th,iz}\) is defined by
Currently the ionization rate \(\langle v_{e} \sigma_{iz}\rangle\) is approximated using the fitting function from [Voronov1997]
where \(A\), \(K\), \(P\) and \(X\) are tabulated for elements up to \(Z=28\). To avoid unphysical negative temperatures, when \(v^2_{th,iz} < 0\) the ionization rate is set to zero in the code.
A similar model of ionization was previously used in Gkeyll and was presented in [Cagas2017].
Charge exchange¶
This process is given by \(i^{+} + n \rightarrow n + i^{+}\), and the simplifed model of this process contained within Gkeyll is based on [Meier2012]. Collision terms appear in the ion and neutral equations as:
where
The cross section is approximated by a fitting function. For hydrogen and Deuterium these are given by, respectively
Wall recycling boundary conditions¶
A model for wall recycling has been implemented at the boundaries where field lines terminate. These boundary conditions provide a source of neutrals from the targets that depends on the flux of outgoing ions. Consider a simulation with one configuration space dimension, parallel to the magnetic field, and three velocity space dimensions (1x3v). Since \(x\) is parallel to the magnetic field, \(v_x\) is the parallel velocity for neutrals. We define the neutral distribution function in the ghost cell at the lower ($x_min$) boundary as
The Maxwellian function for recycled neutrals \(f_{M,rec}\) is defined by a zero mean flow and a temperature that is set in the user input file to model the FranckCondon atoms coming from the wall. The Maxwellian is scaled such that the magnitude of the flux of incoming neutrals is equal to the magnitude of the flux of outgoing ions. At this time, angular dependency is not included in the model of wall recycling.
Neutral species and gyrokinetic plasma species coupling¶
Neutral species are always evolved on the Vlasov grid. For a VlasovMaxwell plasma species, the neutrals and ions are evolved on identical phasespace grids. Thus, the ionneutral interaction terms in Eqs. (2), (4), and (5) are straightforward. However, when the plasma species are evolved using the gyrokinetic model, the ion and neutral velocityspace grids are no longer identical, and it becomes necessary to pass information between two different phasespace grids. This is accomplished by taking fluid moments, \(n\), \(\mathbf{u}\), and \(v^2_{th}\), of the species distribution function and using them to project a Maxwellian distribution function on the destination phasespace grid. This is valid assuming that ion and neutral distribution functions are approximately Maxwellian.
In the VlasovMaxwell formulation, a Maxwellian distribution is defined
where \(d_v\) is the velocityspace dimension. In the gyrokinetic formulation a Maxwellian distribution function is defined
where we have assumed the gyrokinetic grid is either 1X2V or 3X2V. Note that in the gyrokinetic formulation, the fluid velocity moment contains only one component, \(u_\parallel\), which is along the magnetic field line. However, the neutral fluid velocity contains 3 components. It is assumed that once a neutral particle is ionized, the perpendicular components are immediately “smeared out” by the gyromotion. Thus, only the \(z\)component of the neutral fluid velocity moment is included in the Maxwellian projection on the gyrokinetic grid. Conversely, the ion fluid velocity moment contains only one component. Thus, the ion Maxwellian distribution function on the 3V Vlasov grid contains the fluid moment \(\mathbf{u}_i = (u_x = 0, u_y = 0, u_z = u_{\parallel,i})\).
The collision terms in this gyrokineticVlasov coupling become
where \(\mathcal{J}\) is the Jacobian for the gyrokinetic model.
Neutral interactions in Gkeyll input files¶
Electronimpact ionization¶
Below is an example of adding ionization to a VlasovMaxwell simulation:

 App dependencies

local Plasma = (require "App.PlasmaOnCartGrid").VlasovMaxwell()
...
plasmaApp = Plasma.App {

 Common

...

 Species

 VlasovMaxwell electrons
elc = Plasma.Species {
evolve = true,
charge = qe,
mass = me,
...
 Ionization
ionization = Plasma.Ionization {
collideWith = {"neut"},  species to collide with
electrons = "elc",  define name for electron species
neutrals = "neut",  define name for neutral species
elemCharge = eV,  define elementary charge
elcMass = me,  electron mass
plasma = "H",  ion species element
},
...
},
 VlasovMaxwell ions
ion = Plasma.Species {
evolve = true,
charge = qi,
mass = mi,
...
 Ionization
ionization = Plasma.Ionization {
collideWith = {"neut"},  species to collide with
electrons = "elc",  define name for electron species
neutrals = "neut",  define name for neutral species
elemCharge = eV,  define elementary charge
elcMass = me,  electron mass
plasma = "H",  ion species element
},
...
},
 Vlasov neutrals
neut = Plasma.Species {
evolve = true,
charge = 0,
mass = mi,
...
 Ionization
ionization = Plasma.Ionization {
collideWith = {"elc"},  species to collide with
electrons = "elc",  define name for electron species
neutrals = "neut",  define name for neutral species
elemCharge = eV,  define elementary charge
elcMass = me,  electron mass
plasma = "H",  ion species element
},
...
},
},
In order to add ionization to a gyrokinetic simulation and include neutral particles which are evolved using the Vlasov solver, define the Gyrokinetic
App in the dependencies as local Plasma = (require "App.PlasmaOnCartGrid").Gyrokinetic()
. Then replace the neutral Lua table above with
neut = Plasma.Vlasov {
evolve = true,
charge = 0,
mass = mi,
init = Plasma.VmMaxwellianProjection { ... }  initial conditions (and source) defined using Vlasov app
...
 Ionization
ionization = Plasma.Ionization {
collideWith = {"elc"},  species to collide with
electrons = "elc",  define name for electron species
neutrals = "neut",  define name for neutral species
elemCharge = eV,  define elementary charge
elcMass = me,  electron mass
plasma = "H",  ion species element
},
...
bcx = {Vlasov.Species.bcReflect, Vlasov.Species.bcReflect}  boundary conditions defined using Vlasov app
},
Note that Plasma.Species
became Plasma.Vlasov
and Plasma.MaxwellianProjection
became Plasma.VmMaxwellianProjection
but the ionization Lua table remains ionization = Plasma.Ionization
. The latter remains as is since the ionization calculation is carried out from within the Gyrokinetic
App but other parameters such as initial conditions, source, and boundary conditions are defined using the Vlasov
App, which gets called from within the Gyrokinetic
App.
Charge exchange¶
Charge exchange can be added much in the same way as ionization was included above, though the former only affects the ion and neutral species. For the case of gyrokinetic plasma species with Vlasov neutrals, include the following in the Species section of the input file.
 Gyrokinetic ions
ion = Plasma.Species {
evolve = true,
charge = qi,
mass = mi,
...
 Charge exchange
chargeExchange = Plasma.ChargeExchange {
collideWith = {"neut"},  species to collide with
ions = "ion",  define ion species name
neutrals = "neut",  define neutral species name
ionMass = mi,  ion mass
neutMass = mi,  neutral mass
plasma = "H",  ion species element
charge = qi,  species charge
},
...
},
 Vlasov neutrals
neut = Plasma.Vlasov {
evolve = true,
charge = 0,
mass = mi,
...
 Charge exchange
chargeExchange = Plasma.ChargeExchange {
collideWith = {"ion"},  species to collide with
ions = "ion",  define ion species name
neutrals = "neut",  define neutral species name
ionMass = mi,  ion mass
neutMass = mi,  neutral mass
plasma = "H",  ion species element
charge = 0,  species charge
},
...
},
Wall recycling¶
Wall recycling boundary conditions can be included for the neutral Vlasov species by including the following in the neutral table for a simulation in one configurationspace dimension. Since particle fluxes necessary for wall recycling are stored in boundary flux diagnostics, diagnosticBoundaryFluxMoments
must be included in all species table as shown below.
...
 Gyrokinetic electrons
elc = Plasma.Species {
evolve = true,
...
diagnosticBoundaryFluxMoments = {"GkM0", "GkUpar", "GkEnergy"},
}
 Gyrokinetic ions
ion = Plasma.Species {
evolve = true,
...
diagnosticBoundaryFluxMoments = {"GkM0", "GkUpar", "GkEnergy"},
}
 Vlasov neutrals
neut = Plasma.Vlasov {
evolve = true,
charge = 0,
mass = mi,
...
bcx = {Plasma.Vlasov.bcRecycle, Plasma.Vlasov.bcRecycle},
 Recycle elements
recycleTemp = 10*eV,
recycleFrac = 0.5,
recycleIon = "ion",
recycleTime = 100e6,
...
diagnosticBoundaryFluxMoments = {"M0", "u", "M2Flow", "M2Thermal"},
},
Additional flags are required including recycleTemp
which defines \(T_{n,rec}\), recycleFrac
which defines the wall recycling fraction \(\alpha_{rec}\), and recycleIon
which defines the ion species name in the input file. An opptional flag recycleTime
provides a time dependency for the wall recycling fraction, which gradually ramps up to the desired value \(\alpha_{rec,0}\) according to the equation
where \(\tau_{rec}\) is given by recycleTime
.
Examples¶
Two examples of simulations with neutral interactions are presented here. The first uses the VlasovMaxwell solver for the plasma species and includes electronimpact ionization. The second uses the gyrokinetic solver with both electronimpact ionization and charge exchange.
1X1V Vlasov simulation¶
A simple VlasovMaxwell test case in 1X1V with spatially constant fluid moments for all species and periodic boundary conditions can be set up to test conservation properties of this model. Simply run the included input file vlasovIz.lua using standard procedures detailed here. The simulation completes in about 12 seconds on a 2019 MacBook Pro. Then use the Postgkyl commandline tool to check particle and energy conservation. To plot the sum of the integrated particle densities of ions and electrons, use the following command.
pgkyl vlasovIz_ion_intM0.bp vlasovIz_neut_intM0.bp ev 'f[0] f[1] +' plot x 'time' y 'particles'
This produces the plot shown below, illustrating conservation of particle number.
Next plot the sum of integrated thermal energy of ions and neutrals with the following command.
pgkyl vlasovIz_ion_intM2Thermal.bp vlasovIz_neut_intM2Thermal.bp ev 'f[0] f[1] +' plot x 'time' y 'thermal energy'
This produces the plot shown below which demonstrates the conservation of thermal energy.
1X2V gyrokinetic + 1X3V Vlasov simulation¶
This example is based on a simplified model of a scrapeoff layer plasma, the openfield line region in a fusion device. Parameters were chosen based on previous Gkeyll simulations described in [Shi2015]. Gyrokinetic ion and electron species are coupled to Vlasov neutrals via electronimpact ionization and charge exchange interactions. Sheath model boundary conditions are used for the plasma species and reflecting boundary conditions are used for neutrals. The gyrokinetic species are evolved using two velocityspace dimensions, \((v_\parallel, \mu)\). The Vlasov species are run using three velocityspace dimensions, \((v_x, v_y, v_z)\), where the subscripts \((x,y,z)\) correspond to the nonorthogonal fieldline following coordinate system used in the gyrokinetic solver. Thus, \(v_\parallel\) in the gyrokinetic system is identical to the \(v_z\) Vlasov coordinate.
The simulation can be run with the input file 1x2vSOL.lua, which is currently set to run in parallel on 4 processors (decompCuts = {4}
). On a 2019 Macbook Pro, this simulation takes approximately 15 minutes to complete. The output can be analyzed with the Postgkyl tools. For example, the anim
command can be used to observe changes in the electron density profile, as shown below.
pgkyl "1x2vSOL_elc_GkM0_[09]*.bp" interp anim x '$x$' y '$n_e$'
This command produces the following animation of the evolution of the electron density profile in time.
References¶
[Wersal2015]  Wersal, C., & Ricci, P. (2015). A firstprinciples selfconsistent model of plasma turbulence and kinetic neutral dynamics in the tokamak scrapeoff layer. Nuclear Fusion, 55(12), 123014. 
[Voronov1997]  Voronov, G. S. (1997). A practical fit formula for ionization rate coefficients of atoms and ions by electron impact: Z = 128. Atomic Data and Nuclear Data Tables, 65(1), 1–35. 
[Cagas2017]  Cagas, P., Hakim, A., Juno, J., & Srinivasan, B. (2017). Continuum kinetic and multifluid simulations of classical sheaths. Phys. Plasmas, 24(2), 22118. 
[Meier2012]  Meier, E. T., & Shumlak, U. (2012). A general nonlinear fluid model for reacting plasmaneutral mixtures. Physics of Plasmas, 19(7). 
[Shi2015]  Shi, E. L., Hakim, A. H., & Hammett, G. W. (2015). A gyrokinetic onedimensional scrapeoff layer model of an edgelocalized mode heat pulse. Physics of Plasmas, 22(2). 
Projection plugin¶
gkyl Tools¶
Contents
Tools are prepackaged programs that are often useful but not part of the main App system. For example, there are tools to run and query the the regression test system and its output, solve the exact Reimann problem for the Euler equations and compute linear dispersion solvers for multimoment multifluid equations, and coming soon, linear dispersion solvers for kinetic equations.
To get a list of all the supported tools do:
gkyl t
This will give the list of tools and provide a oneline description of what the tool does. Each tool comes with its own embedded help system. For example, to see how to run the exact Euler Reimann solver tool do:
gkyl exacteulerrp h
This will print the help and the command options for the tool to the terminal.
On this page the main tools are documented.
Solving the exact Reimann problem for Euler equations: exacteulerrp
¶
The Reimann problem is a fundamental problem in the theory of
hyperbolic PDEs (like Euler and ideal MHD equations) which brings out
the essential nonlinear structure of the underlying physics. In
particular, the solution to the Reimann problem shows the shocks,
contact discontinuities and rarefactions waves supported by a system
of hyperbolic equations. The exacteulerrp
solves the Reimann
problem exactly for the 1D Euler equations
where
is the total energy contained in the fluid. To get help for this tool do:
gkyl exacteulerrp h
Note that the exacteulerrp
tool itself can print an example input
file to shell that you can then redirect to a file
gkyl exacteulerrp e > sodshock.lua
To run this tool make an input file (or modify the one produced by the e option) with the left and right initial states, the time to compute the solution at and a grid resolution. (The grid is used to determine where the solution is computed. The solution at each grid location is exact and does not depend on the resolution).
lower = 1.0  leftedge of domain
upper = 1.0  rightedge of domain
location = 0.0  location of shock
ncell = 100  Number of cells to use
tEnd = 0.2  Time at which to compute RP
gasGamma = 1.4  Gas adiabatic index
 left/right states { density, velocity, pressure }
leftState = { 1.0, 0.0, 1.0 }
rightState = { 0.125, 0.0, 0.1 }
If this file was saved as “sodshock.lua” you can run the tool:
gkyl exacteulerrp i sodshock.lua
This will produce a set of BP files with the solution to the Riemann problem:
sodshock_density.bp sodshock_pressure.bp
sodshock_internalenergy.bp sodshock_velocity.bp
You can now use postgkyl to plot the solution. For example, to plot all the files in a single figure do:
pgkyl "sodshock_*.bp" pl b
to produce the following plot.
For a comprehensive set of 1D Riemann problems used to benchmark two finitevolume schemes see this note
Linear dispersion relation solver: multimomlinear
¶
The multimomlinear
allows solving linear dispersion equations for
multimoment multifluid equations and will eventually be extended to
full kinetic equations. This tool allows arbitrary number of species,
each of which can be either an isothermal fluid, a fivemoment fluid
or a tenmoment fluid. The fields can be computed from Maxwell
equations or Poisson equations, with the option of some species
“ignoring” the background fields. Certain forms of closures, including
nonlocal HammettPerkins Landau fluid closures, can be used.
For the list of equations and a brief overview of the algorithm used, please see this technical note. Essentially, the key idea of this algorithm is to convert the problem of finding the dispersion relation to an eigenvalue problem and then use a standard linear algebra package (Eigen in this case) to compute the eigensystem. This allows great flexibility as there is no need to directly find complex nonlinear polynomial roots or even formulate the dispersion relation explicitly.
A more careful set of benchmarks and applications are presented in this document.
To run this tool prepare an input file with the species you wish to include, the field and the set of wavenumbers at which the dispersion relation should be computed. An example input file for cold electron and ions is given below.
local Species = require "Tool.LinearSpecies"
 Electrons
elc = Species.Isothermal {
mass = 1.0,  mass
charge = 1.0,  charge
density = 1.0,  number density
velocity = {0.0, 0.0, 0.0},  velocity vector
temperature = 0.0,  temperature
}
 Ions
ion = Species.Isothermal {
mass = 25.0,  mass
charge = 1.0,  charge
density = 1.0,  number density
velocity = {0.0, 0.0, 0.0},  velocity vector
temperature = 0.0,  temperature
}
 EM field
field = Species.Maxwell {
epsilon0 = 1.0, mu0 = 1.0,
electricField = {0.0, 0.0, 0.0},  background electric field
magneticField = {1.0, 0.0, 0.75},  background magnetic field
}
 list of species to include in dispersion relation
speciesList = { elc, ion }
 List of wavevectors for which to compute dispersion relation
kvectors = {}
local kcurr, kmax, NK = 0.0, 4.0, 401
dk = (kmaxkcurr)/(NK1)
for i = 1, NK do
kvectors[i] = {kcurr, 0.0, 0.0}  each kvector is 3D
kcurr = kcurr + dk
end
Any number of species can be specified and the field
can be either
Species.Maxwell
or Species.Poisson
. The wavevectors at which
to compute the dispersion are specified in the kvector
table,
which is a list of three element tables (with components \(k_x,
k_y, k_z\)).
To run this input file (say it is saved in coldplasma.lua):
gkyl multimomlinear i coldplasma.lua
This will create a output file named coldplasma_frequencies.bp, with the eigenvalues stored in a Gkeyll “DynVector” object.
For each element in the dynvector, the first three components are the components of the wavevector and the rest the corresponding \(\omega_n(\mathbf{k})\) with real and imaginary parts stored separately (next to each other). You can plot the real part of the frequencies as function of wavevector (say \(k_x\)) as:
pgkyl coldplasma_frequencies.bp val2coord x0 y 3::2 pl s f0 xlabel "k" ylabel '$\omega_r$' markersize=2
And the imaginary parts as:
pgkyl coldplasma_frequencies.bp val2coord x0 y 4::2 pl s f0 xlabel "k" ylabel '$\omega_r$' markersize=2
Often, it is useful to plot the eigenvalues in the complex plane (real part on Xaxis and imaginary part on the Yaxis). For this do:
pgkyl coldplasma_frequencies.bp val2coord x3::2 y 4::2 pl s f0 xlabel '$\omega_r$' ylabel '$\omega_i$' markersize=2
Note that the frequencies are not outputed in any particular
order. Hence it is not possible to easily extract a single “branch” of
the dispersion relation from the output. Please see pgkyl help to
understand what the val2coord
and pl
(short for plot
) do
and how to use them.
Example of the real frequency for the cold plasma waves is shown below
The species objects can be one of Species.Isothermal
,
Species.Euler
or Species.Tenmoment
. The input parameters
accepted by each of these objects are given below. Note that the input
parameters can either be dimensional or dimensionless. The tool itself
does not use any nondimensionalization.
The Species.Isothermal
takes the following parameters:
Parameter  Description  Default 

mass  Mass of particle  
charge  Charge on particle  
density  Number density  
velocity  Velocity vector {\(v_x\), \(v_y\), \(v_z\)}  
temperature  Fluid temperature (set to zero for cold fluid)  
ignoreBackgroundField  Do not consider background electromagnetic field on species.  false 
The Species.Euler
takes the following parameters:
Parameter  Description  Default 

mass  Mass of particle  
charge  Charge on particle  
density  Number density  
velocity  Velocity vector {\(v_x\), \(v_y\), \(v_z\)}  
pressure  Fluid pressure  
gasGamma  Gas adiabatic index  5/3 
ignoreBackgroundField  Do not consider background electromagnetic field on species.  false 
The Species.Tenmoment
takes the following parameters:
Parameter  Description  Default 

mass  Mass of particle  
charge  Charge on particle  
density  Number density  
velocity  Velocity vector {\(v_x\), \(v_y\), \(v_z\)}  
pressureTensor  Pressure in the fluid frame as a table {\(P_{xx}\), \(P_{xy}\), \(P_{xz}\), \(P_{yy}\), \(P_{yz}\), \(P_{zz}\)}  
useClosure  Flag to turn on various collisionless closures  false 
chi  Multiplicative factor in closure.  \(\sqrt{4/9\pi}\). 
ignoreBackgroundField  Do not consider background electromagnetic field on species.  false 
Note
 Presently, a unmagnetized HammettPerkins closure is implemented. See technote linked above and [Ng2017]. This is not always very useful and accurate, specially for strongly magnetized problems. We hope to implement better closures soon.
 The
ignoreBackgroundField
allows a species to “ignore” the applied background electromagnetic fields. This is often useful when one species is unmagnetized.
References¶
[Ng2017]  J. Ng, A. Hakim, A. Bhattacharjee, A. Stanier, & W. Daughton “Simulations of antiparallel reconnection using a nonlocal heat flux closure”. Physics of Plasmas, 24 (8), 082112–5. (2017) http://doi.org/10.1063/1.4993195 
Postgkyl reference¶
“I love pgkyl! However, Petr, the next feature I want is …” – Ammar Hakim
Below is the living reference to the Postgkyl tool for postprocessing gkyl simulation data. The main aspiration of Postgkyl is to create a postprocessing suit which would:
 Make simple tasks very easy
 Make complex tasks possible
The documentation is split into three sections. The first section introduces the key concepts of using Postgkyl and goes into details of loading data. Particularly the former one is a strongly recommended reading for new users. The second part provides many examples of using Postgkyl, which aim to supplement the examples already presented in the Quickstart chapter. Finally, the last part consists of detailed documentation for each individual available command.
Using Postgkyl¶
Postgkyl is, at its core, a Python library. When properly installed (see Postgkyl install), it can be loaded in any Python script or interactive environment like IPython or JupyterLab.
import postgkyl
Postgkyl can be used to read all the Gkeyll data outputs (including the legacy Gkeyll 1 files), can transform the raw expansion coefficients of Gkeyll basis functions to finitevolume style data, and also contains many postprocessing functions. See the following sections for more details.
Key concepts¶
There are a few basic concepts of using pgkyl
in a command
line. Together, they allow to quickly and easily create quite complex
diagnostics, which would otherwise require writing custom
postpocessing scripts.
Note that there are often multiple ways to achieve the same thing. Sometimes, they are analogous, othertimes, one is superior. This page makes an attempt to explain these key concepts to allow user to choose the best solution for each situation.
Contents
Dataset¶
Each data file which is loaded creates a dataset. Additionally, some commands can create new datasets during the flow.
When files are loaded using wildcard characters, each match creates
its own dataset. Therefore, assuming there are two files, file1.bp
and file2.bp
, located in the current directory, the two following
commands will have the same result; both will create two datasets:
pgkyl file1.bp file2.bp
pgkyl file?.bp
The info command with the c
flag is useful to list
all available datasets.
Command chaining¶
Commands are evaluated from left to right. Each command, by default applies to everything to the left of it. For example, this command chain combines loading data files and interpolate command:
pgkyl file1.bp interpolate file2.bp
file1.bp
is loaded first, its DG expansion coefficients
are then interpolated on a finer uniform mesh, and, finally,
file2.bp
is loaded. interpolate will not be
applied on file2.bp
. This particular example can be used, for
example, to simultaneously work with finiteelement and finitevolume
data.
Commands that should not be applied on all the datasets can be further controlled using tags and by designating some datasets as inactive. Note that there are some commands, e.g., collect, which switch their inputs to inactive themselves.
It is worth noting that there is no limit on how many commands can be chained. See, for example, the particle balance section the the gyrokinetic quickstart page.
Tags¶
The default behavior of most of the commands is agnostic to the tags. For example, the following two commands will lead to the same result:
pgkyl file1.bp file2.bp plot
pgkyl file1.bp t 'f1' file2.bp t 'f2' plot
However, most of the commands can take the use
or u
flag to
limit them only to the datasets with the specified tag. Similar to
the example above, this can be useful when working with different
types of data:
pgkyl file1.bp t 'f1' file2.bp t 'f2' interpolate u f1 plot
Here, interpolate will be used only on the file1.bp
even though it follows loading both of the files. The plot
command
will then apply to both the datasets.
Note that multiple commaseparated tags can be used:
pgkyl file1.bp t 'f1' file2.bp t 'f2' file3.bp t 'f3' interpolate u f1,f2 plot
Additionally, there are some commands like collect or animate are by default tagaware and separate datasets with different tags from each other.
When no tag is specified, the default
tag is assigned.
Warning
When using tags together with wildcard characters, it is important to use quotes, e.g.:
pgkyl 'file?.bp' t name
Without the quotes, the string is replaced with all the matches,
pgkyl
treats them as separate load
commands, and the specified tag is applied only to the last match.
Active and inactive datasets¶
In addition to specifying tags, the flow
of a pgkyl
command chain can be controlled by activating and deactivating
datasets. By default, all loaded datasets are active. This can be
changed with the pair of activate/deactivate and
pg_cmd_deactivate commands. In addition, commands that create a
new dataset, e.g., collect, leave only the output
active. The motivation behind this is that these commands change the
nature of data and user would typically want to keep working only with
the result. The aforementioned collect turns
Ndimensional data to (N+1)dimensional data. With the inputs
inactive, commands can be easily chained, e.g.,
pgkyl 'file*.bp' collect plot
activate/deactivate can either take in indices, tags, or both. When no inputs are specified, everything is activated. The two following commands provide yet another way to to achieve the same result as in the tag example above:
pgkyl file1.bp t f1 file2.bp t f2 activate t f1 interpolate activate plot
pgkyl file1.bp file2.bp activate i 0 interpolate activate plot
In both cases only the file1.bp
is active and, therefore, the
interpolate command is applied only on the first
file. The second activate then reactivates the second file again so
the plot command is going to plot both.
The info command can be useful when working with
multiple active/inactive datasets. Its compact
option shows only
identifiers for each dataset, thus removes some clatter, and
allsets
adds even the currently inactive datasets.
Overwriting vs. new dataset¶
There are two basic ways commandsinteract with inputs. The first type modifies its inputs and pushes data down the chain. A typical example is the interpolate command, which takes expansion coefficients of DG finiteelement data and interpolates them on a finer uniform mesh, essentially creating finitevolume like data.
pgkyl file1.bp interpolate plot
In this case the original information is lost after the interpolate command (lost within this command chain, nothing happens to the data file itself).
The other type does not overwrite its inputs but rather creates a new dataset. As a rule of thumb, these are commands that take (or can take) multiple inputs and/or change the nature of data. Note that these commands often make the result the only active dataset to simplify the flow. A typical example is the ev command:
pgkyl file1.bp file2.bp ev 'f[0] f[1] ' plot
As a result of this chain, there will be three datasets; however, only the result of ev will be active, so the plot command will create just one figure.
There are instances when user does not want to overwrite the
inputs. For example, when we want to use select to
create multiple slices of data. For this purpose, the commands that
would normally overwrite data have the optional tag
or t
flag
which instead creates a new dataset with specified tag. Note that in
this case, the resulting dataset will not be the only one active.
pgkyl file1.bp t input select u input z0 1. t planes \
select u input z0 1. t planes plot u planes
Data loading¶
“I disagree strongly with whatever work this quote is attached to.” – Randall Munroe
One can argue that loading data is the most important part of a
postprocessing tool. In Postgkyl, it is handled by the
postgkyl.data.Data
class (there is a postgkyl.Data
shortcut). It load data on initialization and serves as an input for
all the other parts of Postgkyl.
Examples are provided simultaneously for scripting and command line using output files of an electrostatic twostream instability simulation [twostream.lua].
Contents
Accessing a Gkeyll file¶
Gkeyll files are loaded in Postgkyl by creating a new instance of the
Data
class with the file name as the parameter.
import postgkyl as pg
data = pg.Data('filename')
Next, getGrid()
and getValues()
can be used to return the grid
and values as NumPy arrays. For structured meshes, the getGrid()
return a Python list
of 1D NumPy arrays which represent the nodal
points of the grid in each dimension. Note that since these are nodal
points, these arrays will always have one more cell in each dimension
in comparison to the value array. Another important note is that the
value array always have one extra dimension for
components. Components can represent many things from vector
elements to discontinuous Galerkin expansion coefficients. As a rule,
this extra dimension is always retained even if there is just one
component.
Script example
import postgkyl as pg
data = pg.Data('twostream_elc_0.bp')
print(data.getGrid())
[array([6.283185307179586 , 6.086835766330224 , 5.890486225480862 ,
5.6941366846315 , 5.497787143782138 , 5.301437602932776 ,
5.105088062083414 , 4.908738521234052 , 4.71238898038469 ,
4.516039439535327 , 4.319689898685965 , 4.123340357836604 ,
3.9269908169872414, 3.730641276137879 , 3.5342917352885173,
3.3379421944391554, 3.141592653589793 , 2.945243112740431 ,
2.748893571891069 , 2.552544031041707 , 2.356194490192345 ,
2.1598449493429825, 1.9634954084936211, 1.7671458676442588,
1.5707963267948966, 1.3744467859455343, 1.178097245096172 ,
0.9817477042468106, 0.7853981633974483, 0.589048622548086 ,
0.3926990816987246, 0.1963495408493623, 0. ,
0.1963495408493623, 0.3926990816987246, 0.589048622548086 ,
0.7853981633974483, 0.9817477042468106, 1.178097245096172 ,
1.3744467859455343, 1.5707963267948966, 1.767145867644258 ,
1.9634954084936211, 2.1598449493429825, 2.356194490192344 ,
2.552544031041707 , 2.7488935718910685, 2.9452431127404317,
3.141592653589793 , 3.3379421944391545, 3.5342917352885177,
3.730641276137879 , 3.9269908169872423, 4.123340357836604 ,
4.319689898685965 , 4.516039439535328 , 4.71238898038469 ,
4.908738521234051 , 5.105088062083414 , 5.301437602932776 ,
5.497787143782137 , 5.6941366846315 , 5.890486225480862 ,
6.086835766330225 , 6.283185307179586 ]),
array([6. , 5.8125, 5.625 , 5.4375, 5.25 , 5.0625, 4.875 ,
4.6875, 4.5 , 4.3125, 4.125 , 3.9375, 3.75 , 3.5625,
3.375 , 3.1875, 3. , 2.8125, 2.625 , 2.4375, 2.25 ,
2.0625, 1.875 , 1.6875, 1.5 , 1.3125, 1.125 , 0.9375,
0.75 , 0.5625, 0.375 , 0.1875, 0. , 0.1875, 0.375 ,
0.5625, 0.75 , 0.9375, 1.125 , 1.3125, 1.5 , 1.6875,
1.875 , 2.0625, 2.25 , 2.4375, 2.625 , 2.8125, 3. ,
3.1875, 3.375 , 3.5625, 3.75 , 3.9375, 4.125 , 4.3125,
4.5 , 4.6875, 4.875 , 5.0625, 5.25 , 5.4375, 5.625 ,
5.8125, 6. ])]
print(data.getValues())
[[[ 1.6182154425614533e127 2.2497634664678846e136
2.1705614015952743e127 ... 1.4466223559100639e127
7.7862978418103503e137 2.0112020871650523e136]
[ 7.2163320153412515e118 1.0032681083505769e126
9.6785762877207286e118 ... 6.4497610162539372e118
3.4719259660326997e127 8.9669370964188083e127]
[ 1.3363156717841295e108 1.8578453383418215e117
1.7920360303344134e108 ... 1.1940080895062958e108
6.4284392330301674e118 1.6599988152412963e117]
...
print(data.getGrid()[0].shape)
(65,)
print(data.getGrid()[1].shape)
(65,)
print(data.getValues().shape)
(64, 64, 8)
It is also possible to create an empty instance and fill it using the
push
function.
In the command line mode, a data file is loaded by simply adding it to
the pgkyl
script chain at any position.
pgkyl filename
Note
Under the hood, Postgkyl calls a hidden load
command to load
the file. When provided string does not match any command but is
matching a file, the load command is invoked and the file name is
passed to it. The load command should not be called manually but
it can be used to access the help.
pgkyl load help
Currently, Postgkyl supports h5
file that were used in Gkeyll 1,
Gkeyll 2 ADIOS bp
files, and Gkeyll 0 gkyl
binary files. Many
of the advanced functions like loading only partial data and some
quality of life features like storing the polynomial order of DG
representation are currently available only for the ADIOS bp
files.
Loading multiple datasets¶
Loading multiple files in a script is straightforward; one creates more
instances of the Data
class. Postgkyl does naturally support loading
any number of files.
pgkyl twostream_elc_0.bp twostream_elc_1.bp
All the commands are then generally batch performed on all the data
sets and the plot command creates a separate figure for
each data set (this can be modified with plot options
like f0
).
When batch application of commands is not the desired behavior, some data files can be loaded later in the chain, loaded dataset can be changed from active to inactive (activate/deactivate/pg_cmd_deactivate), or the command scope can be limmited by specifying tags. The Key concepts section provides examples where one desired behavior is achieved in multiple ways. It is left up to the user to chose the preferred one.
Postgkyl also allows for loading with a wild card characters:
pgkyl 'twostream*.bp'
Warning
While the quotes are entirely optional when loading a single file, they change behavior when used with wild card characters. With quotes, a single load command is performed and the wild card matching is done internally by Postgkyl. Without quotes, the wild card is replaced before calling Postgkyl which results in several load command calls. This leads to several key differences:
 With quotes, Postgkyl orders files correctly, i.e.,
file_2
will be beforefile_10
.  With quotes, tags, labels, etc., are applied to all the matching files, not just the last one.
 Some wildcard characters like
[09]
are not supported by every shell.
Using wild card characters might lead to unexpected situations. For
example in the twostream case, the query twostream_elc_*
is
going to return twostream_elc_0.bp
but also the moment files like
twostream_elc_M0_0.bp
. If we want to load just the distribution
functions, we can limit the query. For example:
pgkyl 'twostream_elc_[09]*.bp'
This requires the first character to be a number between 0 and 9, which effectively eliminates all the outputs except for the distribution functions themselves.
Following are details on load parameters which alter the
behavior. Here, we would like to mention that these can be specified
individually for each file of as the global options of the pgkyl
script itself. For example, the partial loading flag z0
(see
bellow) can be applied to one file (file_0
):
pgkyl file_0 z0 0 file_1
Or it can be applied globally to all the files:
pgkyl z0 0 file_0 file_1
This is analogous to:
pgkyl file_0 z0 0 file_1 z0 0
Partial loading¶
Gkeyll output files, especially the higher dimensional ones, can be
large. Therefore, Postgkyl allows to load just a smaller subsection of
each file. This is done with the optional z0
to z5
parameters
for coordinates and comp
for components. Each can be either an
integer number or a string in the form of start:end
. Note that
this does follow the Python convention so the last index is
excluded, i.e., 1:5
will load only the indices/components 1, 2,
3, and 4. This functionality is supported both in the script mode and
the command line mode.
import postgkyl as pg
data = pg.Data('twostream_elc_0.bp', z1='1:3', comp=0)
pgkyl twostream_elc_0.bp z1 1:3 c 0
Note that the select command has a similar use. In addition, it allows to specify a coordinate value instead of an index. However, it requires the whole file to be loaded into memory.
Tags and labels¶
Datasets can be decorated with tags and labels. The former serve mostly to specify the scope of commands (see tags) in the command line mode while the later one allows to add custom labels for plots and printouts.
When no labels are specified, Postgkyl attempts to find the shortest unique identifier and uses it as a label. For example:
pgkyl twostream_elc_0.bp twostream_elc_1.bp info c
0 (default#0)
1 (default#1)
pgkyl twostream_elc_0.bp twostream_field_0.bp info c
elc (default#0)
field (default#1)
pgkyl twostream_elc_0.bp twostream_field_1.bp info c
elc_0 (default#0)
field_1 (default#1)
These labels, can be customized and can include LaTeX syntax, which will be properly rendered in a plot legend.
pgkyl twostream_elc_0.bp l '$t\omega_{pe}=0$' twostream_elc_1.bp l '$t\omega_{pe}=0.5$' info c
$t\omega_{pe}=0$ (default#0)
$t\omega_{pe}=0.5$ (default#1)
Note, the in all these examples, both datasets have the default
tag and are indexed 0
and 1
. These can be manually specified.
pgkyl twostream_elc_0.bp t 'el' twostream_field_0.bp t 'em' info c
elc (el#0)
field (em#0)
Loading data with c2p mapping¶
Warning
This feature was introduced in 1.6.7 and currently only works with
gkyl
binary files.
Postgkyl supports the c2p mapping used in Gkeyll. The file with the
map can be specified using the c2p
keyword. Following are two
plots where a Maxwellian particle distribution is evaluated in
cylindrical coordinates with and without c2p map provided to Postgkyl.
pgkyl rt_eval_on_nodes_fser.gkyl interpolate bms p2 plot a
pgkyl rt_eval_on_nodes_fser.gkyl c2p rt_eval_on_nodes_rthetaten.gkyl interpolate bms p2 plot a
Gkeyll stores the c2p coordinate information as expansion coefficients
of a finite element representation independent of the representation of
the data itself. It is converted to plotting nodal points during the
interpolate command when the information about the data
is provided. However, the interpolate command is never
used when working with finitevolume data. For this instance, the
fv
flag is available which converts the expansion coefficients
to nodal values immediately after loading.
pgkyl euler_axis_sodshockeuler_0.gkyl c2p euler_axis_sodshockmapc2p.gkyl fv select c0 plot a
Note on the colors¶
What makes Postgkyl quite unique, is the wrapping of all the functions into command line commands using the Click Python package. In other words, almost the full functionality of Postgkyl can be used directly in any terminal. This has beneficts for everyday work where typing a single command is faster than writing a Python script and also for work with remote machines and supercomputers, which are primarily accessed through a terminal. However, the classical way of using Postgkyl in a script still provides more control and is well suited, for example, for long and complex scripts for publication level figures.
The command line executable for Postgkyl is called pgkyl
and its
call has the following synopsis:
pgkyl [OPTIONS] COMMAND1 [ARGS]... [COMMAND2 [ARGS]...]...
Here, the COMMAND
represents either data to load or an opperation to be performed on the loaded data. All
the available options can be listed using the inbuilt help, pgkyl
help
or pgkyl h
for short.
Note that Postgkyl supports an arbitrary number of commands. Similar
to piping in Linux, results of one command are passed to another. This
allows Postgkyl to perform even a rather complex diagnostics and
create complicated plots straight up from the terminal. However, this
puts a responsibility on the user to ensure that the commands are
called in a logical order. Similarly to the main part of pgkyl
,
help
can be called for each individual command which will
provide additional information. Finally, it is worth mentioning that
it is not necessary to write the full names of each command and the
shortest unique sequence is good enough. Still, full names will be
used through this documentation for clarity.
Examples¶
The Gkeyll workflow typically involves a mixture of command line and script mode use of postgkyl. It is useful to become versed in both of these, and to that end we provide a number of examples below.
Commandline examples¶
There are more examples in the pages describing each of the pgkyl
commands and functions. For a brief set of examples we here consider
the output of a twostream instability VlasovMaxwell simulation.
List outputs¶
We can list the different kinds of files outputted by the simulation with the listoutputs command
pgkyl list
twostream_elc
twostream_elc_M0
twostream_elc_M1i
twostream_elc_M2ij
twostream_elc_M3i
twostream_field
Notice that we used the abbreviation list
instead of listoutputs
;
this is allowed and whenever convenient we will use abbreviations if
the meaning is clear. This writes out a list of the unique filename
stems of the files outputted by a simulation.
Obtain file information¶
Now suppose we wish to know more about one file in particular, say the electron initial condition. We can load the data set and probe it with the info command:
pgkyl twostream_elc_0.bp info
Set (default#0)
├─ Time: 0.000000e+00
├─ Frame: 0
├─ Number of components: 8
├─ Number of dimensions: 2
├─ Grid: (uniform)
│ ├─ Dim 0: Num. cells: 64; Lower: 6.283185e+00; Upper: 6.283185e+00
│ └─ Dim 1: Num. cells: 64; Lower: 6.000000e+00; Upper: 6.000000e+00
├─ Maximum: 3.804653e+00 at (31, 26) component 0
├─ Minimum: 6.239745e01 at (31, 38) component 2
├─ DG info:
│ ├─ Polynomial Order: 2
│ └─ Basis Type: serendipity (modal)
├─ Created with Gkeyll:
│ ├─ Changeset: 82231b06678c+
│ └─ Build Date: Feb 7 2021 09:14:06
The output shows (in order):
 Simulation time of this snapshot.
 Frame number of this snapshot.
 Number of degrees of freedom per cell (components). In this case 8 basis functions for the second order 2D Serendipity basis ([Arnold2011], see Modal basis functions for more details).
 Number of dimensions (in this case 2 because it is a 1x1v simulation, 1D in position and 1D in velocity).
 The grid resolution and extents.
 The maximum value of the dataset (in this case largest DG coefficient).
 The minimum value of the dataset (in this case smallest DG coefficient).
 The polynomial order and type of basis.
 The Gkeyll build used to produce this simulation.
Plotting (interpolated data or single coefficients)¶
In the case of DG data one is most frequently interested in plotting
the actual function instead of its expansion (DG) coefficients. To do so
we construct finitevolumelike dataset on a uniform mesh with higher
resolution using the interpolate command. If we interpolate
the same dataset (electron initial condition) as in the previous example
and inspect it with info
, we get
pgkyl twostream_elc_0.bp interpolate info
Set (default#0)
├─ Time: 0.000000e+00
├─ Frame: 0
├─ Number of components: 1
├─ Number of dimensions: 2
├─ Grid: (uniform)
│ ├─ Dim 0: Num. cells: 192; Lower: 6.283185e+00; Upper: 6.283185e+00
│ └─ Dim 1: Num. cells: 192; Lower: 6.000000e+00; Upper: 6.000000e+00
├─ Maximum: 1.970512e+00 at (96, 79)
├─ Minimum: 1.177877e08 at (94, 60)
├─ DG info:
│ ├─ Polynomial Order: 2
│ └─ Basis Type: serendipity (modal)
├─ Created with Gkeyll:
│ ├─ Changeset: 82231b06678c+
│ └─ Build Date: Feb 7 2021 09:14:06
Note that now there is only one degree of freedom (component) per cell but there are 3X more cells; we have interpolated the function onto a finer mesh.
We can then plot the distribution function evaluated on this finer mesh with the plot command
pgkyl twostream_elc_0.bp interp plot
Note the allowed abbreviation of interpolate
to interp
. This command
produces the figure below:
In some cases one may also be interested in plotting single expansion
coefficients. The most common use is to, for example, plot the cellaveraged
value (which is the zeroth coefficient times a constant). We can select
specific coefficients in all cells with the select and the
c
flag. So, to plot the cell averaged electron distribution function
(times a factor) we would use
pgkyl twostream_elc_0.bp select c0 plot
Which produces the figure below. Notice how the values are slightly different and the resolution is coarser than in the previous plot of interpolated data.
Plot data slices¶
The select command introduced in the previous example can also be
used with the z
flag in order to select a data slice which we may
subsequently plot. In the two stream instability example, the electron initial
condition is clearly independent of position, so we can plot as a function of
velocity space at \(x=2.0\) with
pgkyl twostream_elc_0.bp interp select z0 2. plot
producing the following plot
Postgkyl is currently limited to 1D and 2D plots, so in order to visualize
datasets that have more than 3 dimensions one may need to select several slices
at once. You can do that with multiple z
flags. For example, if we have a
2x2v simulation producing fourdimensional distribution functions, we can select
a slice at \(v_x=0.\) and \(v_y=1.\) with z2 0. z3 1.
.
Create animations¶
Another useful operation is to load multiple datasets consisting of consecutive
frames in a timedependent simulation, plot them and stitch them together into a
movie. This can be accomplished by the animate command. We first
load all the frames containing the electron momentum densities in time
(i.e. elc_M1i), then interpolate them onto a finer mesh, and put together all the
frames with the animate
command:
pgkyl "twostream_elc_M1i_[09]*.bp" interp animate x 'x' y 'Momentum'
Here we are also using the x
and y
flags to the animate
command in
order to place labels in the figure. The product is the movie given below:
Of course, one can use command chaining to slice higher dimensional data prior to
calling the animate
command. For example, creating an animation of the
distribution function along velocityspace at \(x=0\) would be accomplished with
pgkyl "twostream_elc_[09]*.bp" interp sel z0 0. animate
Here we have used the abbreviations interp
and sel
in favor of
interpolate
and select
, respectively. Such command produces this animation:
Averaging and integrating¶
A common diagnostic need is to perform averages and integrations over a
dimension or over time. In general averages can be performed either by
using the ev command with the avg
operation, or using
the integrate and later dividing by the corresponding
space/time segment (with the ev command).
We use the final electron distribution function (frame 100) as an example. Let’s first plot it in phasespace to get a sense of it:
pgkyl twostream_elc_100.bp interp pl
Suppose we wish to average over the central region \(x\in[2,2]\) where a velocityspace ‘hole’ forms. We can plot such \(x\) integral as follows
pgkyl twostream_elc_100.bp interp sel z0 2.:2. ev 'f[0] 0 avg' pl
As mentioned above, we can also do this with the integrate command. We accomplish that with
pgkyl twostream_elc_100.bp interp sel z0 2.:2. integrate 0 ev 'f[0] 4. /' pl
Another useful application of ev
(with avg
) and integrate
is to average or integrate quantities over time. Consider the evolution of
the electron distribution function along velocity space at \(x=0\) in
the previous example. The action starts after the 60th frame approximately,
so if we wish to timeaverage the distribution function
at \(x=0\) we could use the ev command as follows:
pgkyl "twostream_elc_[09]*.bp" activate i '59:' interp sel z0 0. collect ev 'f[0] 0 avg' pl
A different way to accomplish the same time average over frames 59100 and dividing by the corresponding time period (\(\tau=5029.5018=20.4982\)):
pgkyl "twostream_elc_[09]*.bp" activate i '59:' interp sel z0 0. collect integrate 0 ev 'f[0] 20.4982 /' pl
To break this last approach down, the command does the following (in order):
 Load all frames of the electron distribution function.
 Activate only frames 59100.
 Interpolate each frame onto a finer mesh and slice at \(x=0\).
 Collect all the slices into a single dataset. This produces a 2D dataset with time along the 0th dimension and velocityspace along the 1st dimension.
 Use the integrate command to integrate along the 0th dimension (time).
 Use ev to divide the timeintegrated quantity by the appropriate time window.
 Plot.
The product of either of these comands is shown below:
Plot differences between datasets¶
It is common to have to evaluate the difference between two datasets. These could be two frames from the same simulation, or two datasets from different simulations. There are also various ways to discern differences, and below we show how to plot them in a single figure or how to plot their actual difference.
Suppose we wish to see how the electron distribution function has changed along \(x\) between \(t=0\) (0th frame) and \(t=50\) (100th frame) at \(v_x=0\). We can plot both of these datasets as follows
pgkyl twostream_elc_0.bp twostream_elc_100.bp interp sel z1 0. plot f0 x '$x$' y '$f_e(x,v_x=0)$'
where we have specified the figure with f0
so they are both plotted
together, and we have used x
and y
to place labels. The plot
produced is
Another alternative is to compute the actual difference of the two data sets with ev:
pgkyl twostream_elc_0.bp twostream_elc_100.bp interp sel z1 0. ev 'f[1] f[0] ' pl
Note that these operations also work with 2D datasets. So we could’ve have taken the whole distribution function in phase space at \(t=0,50\), subtracted them and plot them with
pgkyl twostream_elc_0.bp twostream_elc_100.bp interp ev 'f[1] f[0] ' plot diverging
which thanks to the diverging
flag, produces the following image:
Saving plots/animations to a file¶
Any of the figures above can be saved to a file by appending either save
,
or saveas
followed by the desired filename. For example the diverging
2D colorplot in the previous section can be saved to a file with
pgkyl twostream_elc_0.bp twostream_elc_100.bp interp ev 'f[1] f[0] ' plot diverging saveas 'twostream_elc_fr0m100.png'
Fileformats supported depend on matplotlib, but likely include .png, .pdf and .eps.
The animate command is also capable of saving animations to a file. It requires an ffmpeg installation, and once that is available we can create an animation described earlier as
pgkyl "twostream_elc_[09]*.bp" interp sel z0 0. animate saveas 'twostream_elc_z0eq0p0.mp4'
Extracting input file¶
In our commitment to reproducibility, Gkeyll output files store the Lua input file used to produce that data. This input file can be extracted using the extractinput command, as follows
pgkyl twostream_elc_0.bp extractinput
By default, this commands simply prints the input file to screen. However, this could be easily piped into a new file with
pgkyl twostream_elc_0.bp extractinput > newInputFile.lua
Reference¶
[Arnold2011]  Arnold, D. N. and Awanou, G. “The serendipity family of finite elements.” Foundations of Computational Mathematics 11.3 (2011): 337344. 
Script mode examples¶
Examples of using the script mode of postgkyl.
Simple loading of data in grids:
import postgkyl as pg
def func_data(ionDensityData):
ionDensityInterp = pg.data.GInterpModal(ionDensityData, 1, 'ms')
interpGrid, ionDensityValues = ionDensityInterp.interpolate()
# get cell center coordinates
CCC = []
for j in range(0,len(interpGrid)):
CCC.append((interpGrid[j][1:] + interpGrid[j][:1])/2)
x_vals = CCC[0]
y_vals = CCC[1]
z_vals = CCC[2]
X, Y = np.meshgrid(x_vals, y_vals)
ionDensityGrid = np.transpose(ionDensityValues[:,:,z_slice,0])
return x_vals,y_vals,X,Y,ionDensityGrid
ionDensity=DIR+'<filename>_ion_GkM0_%d'%data_num[i]+'.bp'
ionDensityData = pg.data.GData(ionDensity)
x_vals,y_vals,X,Y,ionDensityGrid = func_data(ionDensityData)
Evaluating derivative of data on grids:
import postgkyl as pg
def func_data(phiData):
phiInterp = pg.data.GInterpModal(phiData, 1, 'ms')
interpGrid, phiValues = phiInterp.interpolate()
exValues =  np.gradient(phiValues,dx,axis = 0)
dexdxValues = np.gradient(exValues,dx,axis = 0)
eyValues =  np.gradient(phiValues,dy,axis = 1)
# get cell center coordinates
CCC = []
for j in range(0,len(interpGrid)):
CCC.append((interpGrid[j][1:] + interpGrid[j][:1])/2)
x_vals = CCC[0]
y_vals = CCC[1]
z_vals = CCC[2]
X, Y = np.meshgrid(x_vals, y_vals)
exGrid = np.transpose(exValues[:,:,z_slice,0])
eyGrid = np.transpose(eyValues[:,:,z_slice,0])
dexdxGrid = np.transpose(dexdxValues[:,:,z_slice,0])
return x_vals,y_vals,X,Y,exGrid,eyGrid,dexdxGrid
phi=DIR+'<filename>_phi_%d'%data_num[i]+'.bp'
phiData = pg.data.GData(phi)
x_vals,y_vals,X,Y,exGrid,eyGrid,dexdxGrid = func_data(phiData)
Function/Command reference¶
This is the reference for Postgkyl functions and commands. In the most cases, the commands are simple click wrappers of the Python function. However, where a command addresses something that is easilly done in a script, there is no Python function equivalent.
Examples are provided simultaneously for the analogous functions and commands using output files of an electrostatic twostream instability simulation [twostream.lua].
animate¶
Create movies (animations) by stitching together figures created by plotting multiple datasets, typically originating from a data load with wildcard/regex.
Saving the animation to a file requires an ffmpeg installation.
Command line¶
Command help
pgkyl animate h
Usage: pgkyl anim [OPTIONS]
Animate the actively loaded dataset and show resulting plots in a loop.
Typically, the datasets are loaded using wildcard/regex feature of the f
option to the main pgkyl executable. To save the animation ffmpeg needs to
be installed.
Options:
u, use TEXT Specify a tag to plot.
s, squeeze Squeeze the components into one panel.
b, subplots Make subplots from multiple datasets.
nsubplotrow INTEGER Manually set the number of rows for subplots.
nsubplotcol INTEGER Manually set the number of columns for subplots.
transpose Transpose axes.
c, contour Make contour plot.
q, quiver Make quiver plot.
l, streamline Make streamline plot.
g, group [01] Switch to group mode.
s, scatter Make scatter plot.
markersize FLOAT Set marker size for scatter plots.
style TEXT Specify Matplotlib style file (default:
Postgkyl).
d, diverging Switch to diverging colormesh mode.
arg TEXT Additional plotting arguments, e.g., '*'.
a, fixaspect Enforce the same scaling on both axes.
logx Set xaxis to log scale.
logy Set yaxis to log scale.
logz Set values of 2D plot to log scale.
xscale FLOAT Value to scale the xaxis (default: 1.0).
yscale FLOAT Value to scale the yaxis (default: 1.0).
vmax FLOAT Set maximal value of data for plots.
vmin FLOAT Set minimal value of data for plots.
f, float Choose min/max levels based on current frame
(i.e., each frame uses a different color range).
xlim TEXT Set limits for the xcoordinate (lower,upper)
ylim TEXT Set limits for the ycoordinate (lower,upper).
legend / nolegend Show legend.
forcelegend Force legend even when plotting a single
dataset.
x, xlabel TEXT Specify a xaxis label.
y, ylabel TEXT Specify a yaxis label.
clabel TEXT Specify a label for colorbar.
title TEXT Specify a title.
i, interval INTEGER Specify the animation interval.
save Save figure as PNG.
saveas TEXT Name to save the plot as.
fps INTEGER Specify frames per second for saving.
dpi INTEGER DPI (resolution) for output.
e, edgecolors TEXT Set color for cell edges (default: None)
showgrid / noshowgrid Show gridlines (default: True)
collected Animate a dataset that has been collected, i.e.
a single dataset with time taken to be the first
index.
hashtag Turns on the pgkyl hashtag!
show / noshow Turn showing of the plot ON and OFF (default:
ON).
h, help Show this message and exit.
Suppose we extend the simulation time of the
gyrokinetic ion acoustic wave simulation
to tEnd=50
and the number of frames to nFrame=100
, and run it
with
gkyl gkionSound1x2vp1.lua
Note
This simulation took 582 seconds running with decompCuts=8
(8 cores)
on a 2020 MacBookPro.
and plot the electrostatic potential from all frames with
pgkyl "gkionSound1x2vp1_phi_[09]*.bp" interp collect pl x '$x$' y '$\phi$' group 1 clabel 'time'
this produces the pretty picture below, showing the electrostatic potential as a function of \(x\) in each frame, with the color indicating the time stamp of that frame.
This kind of plot can become impractical if there are many frames, or if one would instead like to see a timedependent movie of the evolution of the potential. For that purpose we can create an animation of \(\phi(x,t)\) using the following command
pgkyl "gkionSound1x2vp1_phi_[09]*.bp" interp anim x '$x$' y '$\phi$' saveas 'gkionSound1x2vp1_phi.mp4'
Note
In order to save the animation to an .mp4 file, simply append
saveas fileName.mp4
to the end of the above command.
and this will produce the animation below:
One can clearly see the amplitude of the wave decaying as Landau damping ruins the fun.
collect¶
Assemble multiple active datasets into a new, single dataset.
It is also possible to collect datasets into chunks (multiple datasets) of a specified size rather than into a single dataset.
Command line¶
pgkyl collect help
Usage: pgkyl collect [OPTIONS]
Collect data from the active datasets and create a new combined dataset.
The timestamp in each of the active datasets is collected and used as the
new Xaxis. Data can be collected in chunks, in which case several
datasets are created, each with the chunksized pieces collected into each
new dataset.
Options:
s, sumdata Sum data in the collected datasets (retain components)
p, period FLOAT Specify a period to create epoch data instead of time
data
offset FLOAT Specify an offset to create epoch data instead of time
data [default: 0.0]
c, chunk INTEGER Collect into chunks with specified length rather than
into a single dataset
u, use TEXT Specify a 'tag' to apply to (default all tags).
t, tag TEXT Specify a 'tag' for the result.
l, label TEXT Specify the custom label for the result.
h, help Show this message and exit.
There are many uses of the collect
command. One such use, for
example, is to plot quantities over time. If we consider the
two stream instability VlasovMaxwell simulation
we can gather the \(x\)component of the electric field from
all frames and plot it as a function of time and space with
pgkyl "twostream_field_[09]*.bp" interp sel c0 collect pl x 'time' y 'x'
The field is very small at first but is exponentially growing, so in the above figure we mostly see the saturated electric field at later times.
The collect command is also able to collect datasets into chunks
of a specified length with the c
flag. So suppose we wish to
compute a time average of the electric field over three separate
periods in the two stream instability simulation (we will set
nFrame=98
in the gkyl input file for this example so there are
a total of 99 frames), we can collect the field data into three
separate datasets with 33 frames each and use the ev
command to average in time (0th dimension) as follows:
pgkyl "twostream_field_[09]*.bp" interp sel c0 collect c 33 ev l 'frames 032' 'f[99] 0 avg' \
ev l 'frames 3365' 'f[100] 0 avg' ev l 'frames 6698' 'f[101] 0 avg' activ i102:105 \
pl x '$x$' y '$\langle E_x\rangle_t$'
producing the following three plots
Clearly the field amplitude averaged over 30 frame intervals is increasing as the instability develops.
Finally, collect
allows a transformation of the time dimension
so that instead of time t
it instead becomes (toffset)/period
via the offset
and period flags
. This is used,
for example, in astronomy for variable stars and in creating
Poincare plots (see http://ammarhakim.org/sj/je/je32/je32vlasovtestptcl.html).
current¶
Command line¶
Command help
pgkyl velocity h
Usage: pgkyl curr [OPTIONS]
Accumulate current, sum over species of charge multiplied by flow
Options:
q, qbym TEXT Flag for multiplying by charge/mass ratio instead of just
charge [default: False]
u, use TEXT Specify a 'tag' to apply to (default all tags).
t, tag TEXT Tag for the resulting current array [default: current]
l, label TEXT Custom label for the result [default: J]
h, help Show this message and exit.
differentiate¶
Compute the derivatives of a function along each direction.
This returns a dataset with as many components as there are
dimensions in the original dataset, i.e. it returns the gradient,
unless the d
flag is used to specify a direction.
This command takes the basis expansion in a cell, differentiates it analytically in the desired directions, and then interpolates it onto a finer mesh in the same fashion as the interpolate command.
Note that differentiation is also possible with the ev
command and the grad
operation.
Command line¶
Command help
pgkyl differentiate h
Usage: pgkyl differentiate [OPTIONS]
Interpolate a derivative of DG data on a uniform mesh
Options:
b, basistype [msnsmo] Specify DG basis
p, polyorder INTEGER Specify polynomial order
i, interp INTEGER Interpolation onto a general mesh of specified
amount
d, direction INTEGER Direction of the derivative (default: calculate
all)
r, read BOOLEAN Read from general interpolation file
t, tag TEXT Specify a 'tag' to apply to (default all tags).
o, outtag TEXT Optional tag for the resulting array
h, help Show this message and exit.
Let’s take a gyrokinetic ion acoustic wave simulation as an example. We can examine the initial electrostatic potential generated by the initial conditions with
pgkyl gkionSound1x2vp1_phi_0.bp interp pl x '$x$' y '$\phi$'
giving the plot shown below on the left.
Suppose we wished to know what the initial electric field is, then we would differentiate the potential and multiply it by 1 as follows
pgkyl gkionSound1x2vp1_phi_0.bp diff d 0 ev 'f[1] 1 *' pl x '$x$' y '$\partial_x\phi$'
Note we we have abbreviated differentiate
> diff
, either is allowed.
This command produces the electric field above on the right. It is cellwise
constant because we use a piecewise linear basis function.
Now suppose you wish to examine the gradient of the ion distribution function at \(t=0\) and \(x=0\). This can be accomplished with the following command
pgkyl gkionSound1x2vp1_ion_0.bp diff sel z0 0. pl x '$v_\parallel$' y '$\mu$'
in order to produce:
Starting with the top left and going clockwise, this plot provides the gradient in \(f_i(x,v_\parallel,\mu)\) along \(x\), \(v_\parallel\) and \(\mu\), all three evaluated at \(x=0\).
euler¶
Command line¶
Command help
pgkyl euler help
Usage: pgkyl euler [OPTIONS]
Compute Euler (fivemoment) primitive and some derived variables from
fluid conserved variables.
Options:
u, use TEXT Specify a 'tag' to apply to (default all
tags).
g, gas_gamma FLOAT Gas adiabatic constant
v, variable_name [densityxvelyvelzvelvelpressurekemach]
Variable to extract
h, help Show this message and exit.
Script mode¶
ev¶
ev
is a special command that enables simple operations like adding
or multiplying datasets directly in a terminal. As these operations
are generally simple to do in a script mode, ev
is available only
in the command line interface.
Contents
Reverse Polish notation¶
ev
uses reverse Polish notation (RPN; also know as postfix
notation) to input operations. Unlike in the more common infix
notation, in RPN, the operators follow operands. For example 4  2
is written in RPN as 4 2 
; 4  (2 / 5)
is 4 2 5 / 
. This
can be further demonstrated on the sequence that processes this
expression.
4
>>>4
(Add 4 to the stack)4 2
>>>4 2
(Add 2 to the stack)4 2 5
>>>4 2 5
(Add 5 to the stack)4 2 5 /
>>>4 0.4
(Divide the last two elements with each other)4 2 5 / 
>>>3.6
(Subtract the last two elements from each other)
More information about RPN can be found, for example, on Wikipedia.
Using ev
on simple datasets¶
The ev
command has one required argument, which is the operation
sequence. Datasets are specified in the sequence with
f[set_index][component]
. Note that Python indexing conventions are
used. Specifically, when no indices are specified, everything is used,
i.e., f[0][:]
and f[0]
are treated as identical calls. The
indices need to be integers unless a global indexing mode g
is
used, which ignores active and inactive sets. Then, Python
slice syntax can be used as well including negative indices and
strides. For example, f[2:1:2]
will select every other dataset,
starting with the third one (zeroindexed), and ending with one before
the last (by Python conventions, the upper bound is not included; in
this example it might end with the third element from the end because
of the stride).
Note
Dataset selection in ev
internally uses the same code as the
activate/deactivate/pg_cmd_deactivate commands, so the
following commands produce similar results (ev
actually copies
the datasets instead of just activating some).
pgkyl twostream_elc_?.bp activate i '2:1:2'
pgkyl twostream_elc_?.bp ev g 'f[2:1:2]'
However, this is probably a fringe application of ev
.
The simplest example of ev
is a numerical operation performed on
a dataset, e.g., dividing the values by the insidious factor of 2:
pgkyl twostream_elc_0.bp ev 'f[0] 2 /'
This can be also combined with the fact that ev
can access dataset
metadata as long as they are included (which is a new feature in
Gkeyll introduced in January 2021). An example of this can be plotting
number density from a fluid simulation (Gkeyll outputs mass density).
pgkyl 5m_fluid_elc_0.bp ev 'f[0][0] f[0].mass /' plot
Note that on top of dividing by mass, only the first component, which corresponds to density, was selected. This can be easily extended to apply on multiple datasets and create an animation using the animate command
pgkyl '5m_fluid_elc_[09]*.bp' ev g 'f[:][0] f[:].mass /' animate
The capabilities are not limited to operations with float factors. As
an example, ev
can be used to visualize differences
(diverging
mode of the plot command is well suited
for this)
pgkyl twostream_elc_0.bp twostream_elc_80.bp interpolate ev 'f[1] f[0] ' plot diverging
Note
info command, especially with the compact
c
flag can be useful to print indices for available datasets.
The same concept can be used to calculate bulk velocity from the first two moments:
pgkyl twostream_elc_M0_0.bp twostream_elc_M1i_0.bp interpolate ev 'f[1] f[0] /' plot
Finally, it is worth noting that this syntax cannot be used when there are datasets with more than one tag active.
Using ev
on datasets with tags¶
The ev
command is tagaware. Tagged datasets use the following
notation t.tag_name[set_index][component]
. Using this, the
previous example can be reproduced:
pgkyl twostream_elc_M0_0.bp t dens twostream_elc_M1i_0.bp t mom interp ev 'mom[0] dens[0] /' plot
However, unlike the previous example, this can be naturally extended for batch loading and animate:
pgkyl 'twostream_elc_M0_[09].bp' t dens 'twostream_elc_M1i_[09]*.bp' t mom interp ev g 'mom dens /' animate
Examples of specific ev
operations¶
In this section we provide examples of some ev
operations that
are less trivial or intuitive.
grad¶
This operation differentiates a along a direction given by the second operand. So, for example, given the data from an ion sound wave gyrokinetic simulation we can plot the initial electrostatic potential with
pgkyl gkionSound1x2vp1_phi_0.bp interp pl x '$x$' y '$\phi$'
and compute the parallel electric field by differentiating the potential along \(x\) as follows:
pgkyl gkionSound1x2vp1_phi_0.bp interp ev 'f[0] 0 grad 1 *' pl x '$x$' y '$\phi$'
These produce the following plots:
int¶
Integrate a dataset along a direction specified by the second operand,
or along multiple directions specified by a commaseparated list. If we
once again take the
ion sound wave gyrokinetic simulation
data, we can examine the number of particles in the simulation (should be
conserved) by taking the time trace of the integrated ion number density
(intM0
) and taking its mean:
pgkyl gkionSound1x2vp1_ion_intM0.bp ev 'f[0] mean' pr
which prints out
12.566370614358858
If we instead use ev
to integrate the initial and/or the final number
density GkM0
, we should get roughly the same answer. We can check that
this is the case by typing
pgkyl gkionSound1x2vp1_ion_GkM0_10.bp interp ev 'f[0] 0 int' pr
which produces
12.566370614358522
and we have shown that the number of particles at the end is roughly the same as the mean number of particles throughout the simulation.
extractinput¶
Data files produced by Gkeyll may have the input file used
to generate that data saved in them as an attribute. We can
extract the input file and print it to screen or to a new
file with extractinput
.
Certain versions of ADIOS (the I/O library Gkeyll uses) have
trouble saving very long files. For this reason saving the
input file in data files is sometimes disabled. In those cases
extractinput
may fail or return nothing.
Command line¶
Command help
pgkyl extractinput h
Usage: pgkyl extractinput [OPTIONS]
Extract embedded input file from compatible BP files
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
h, help Show this message and exit.
If we run the ion acoustic gyrokinetic simulation,
it would produce save the initial ion distribution function in
the file gkionSound1x2vp1_ion_0.bp
. We can extract the
input file used to generate this data with
pgkyl gkionSound1x2vp1_ion_0.bp extract
where we have shortened extractinputfile
into extract
for simplicity. This would print it to screen.
Suppose you wish to write it to a new file because, say, you lost the original input file and wish to recreate or restart this simulation. You could write the saved input file to a new input file with:
pgkyl gkionSound1x2vp1_ion_0.bp extract > gkionSound1x2vp1.lua
If you are using the same name as the original input file you’ll be able to restart the simulation. If the original input file is still in the directory it will be overwritten.
fft¶
Command line¶
Command help
pgkyl fft help
Usage: pgkyl fft [OPTIONS]
Calculate the Fourier Transform or the powerspectral density of input
data. Only works on 1D data at present.
Options:
p, psd Limits output to positive frequencies and returns the power
spectral density FT^2.
i, iso Bins power spectral density FT^2, making 1D power
spectra from multidimensional data.
u, use TEXT Specify a 'tag' to apply to (default all tags).
t, tag TEXT Optional tag for the resulting array
l, label TEXT Custom label for the result
h, help Show this message and exit.
Fourier analysis on onedimensional data is available in pgkyl. This does not mean that the simulation has to be 1D, but for higherD simulations one must first reduce the dimensionaly of the data before doing an FFT.
For example, we take a simple incompressible Euler 2D simulation of a velocity shear instability (Kelvin Helmholtz). Using the pgkyl command
pgkyl "incompEulerKH2xp1_fluid_[09]*.bp" interp anim a x 'x' y 'y' clabel '$\psi(x,y)$'
we obtain the following velocity potential
This simulation was initialized with a slow variation in x and a small but more oscillatory perturbation in y:
where \(\alpha=0.05\), \(k_x=2\pi/L_x\), \(k_y=16\pi/L_y\), and \(L_y=4L_x=40\). If we were to examine the Fourier transform of this initial condition at \(x=5\) with
pgkyl incompEulerKH2xp1_fluid_0.bp interp sel z0 5. fft ev 'f[0] abs' pl x '$k_y$' y '$\psi_{k_y}(x=5)$' logy xscale 6.283185
where we scaled the \(x\)axis by \(2\pi\) because of SciPy’s fftfreq convention, we would obtain
or most commonly one looks at the power spectrum of a signal, which we can obtain
with the p
flag:
pgkyl incompEulerKH2xp1_fluid_0.bp interp sel z0 5. fft p pl x '$k_y$' y '$\psi_{k_y}(x=5)^2$' xscale 6.283185 logy
This plot has the peak we would expect at \(k_y=16\pi/L_y=1.2566\), but it also has two other peaks we did not expect. This is because we are FFTing interpolated DG data which introduces modes if the transform is not done weakly (not covered here).
We could also look at how this spectrum changes in time with the following command
pgkyl "incompEulerKH2xp1_fluid_[09]*.bp" interp sel z0 5. fft p collect pl group 1 logy x '$k_y/(2\pi)$' y '$\left\psi_{k_y}(x=5)\right^2$' clabel 'time'
showing how the spectrum goes from being peaked at specific \(k_y\)’s to being a fully filled spectrum when turbulence sets in.
Script mode¶
Although it is possible to call postgkyl’s fft command from a script, we recommend that you instead use SciPy’s fft package directly (or alternatively NumPy’s fft package). Here’s an example of how to perform the same FFT described in the command line section above
import postgkyl as pg
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt
import numpy as np
fileName = 'incompEulerKH2xp1_fluid_0.bp'
polyOrder, basisType = 1, 'ms'
pgData = pg.GData(fileName)
pgInterp = pg.GInterpModal(pgData, polyOrder, basisType)
pgInterp.interpolate(overwrite=True)
yInt, psi_z0eq5p0 = pg.data.select(pgData, z0=5.0)
y = 0.5*(yInt[1][:1]+yInt[1][1:])
ky, psi_z0eq5p0_ky = 2.*np.pi*fftfreq(np.size(y)), fft(np.squeeze(psi_z0eq5p0))
plt.semilogy(ky, np.abs(psi_z0eq5p0_ky))
plt.show()
where we used the select command to pick the data at \(x=5\), transformed the \(y\) coordinates from nodal to cellcenter coordinates, and squeezed the data to remove redundant dimensions.
growth¶
Measure the growth rate in the time trace of a quantity, assuming that it is a squared quantity (e.g. \(\phi^2\), \(\mathbf{E}^2\)). If one applies this to a nonsquared quantity then the growthrate would be off by a factor of two.
The command works by fitting \(A e^(2\gamma t)\) to the data and returning \(\gamma\).
Command line¶
Command help
$ pgkyl growth help
Usage: pgkyl growth [OPTIONS]
Attempts to compute growth rate (i.e. fit e^(2x)) from DynVector data,
typically an integrated quantity like electric or magnetic field energy.
Options:
u, use TEXT Specify a 'tag' to apply to (default all
tags).
g, guess <FLOAT FLOAT>... Specify initial guess
p, plot Plot the data and fit
minn INTEGER Set minimal number of points to fit
maxn INTEGER Set maximal number of points to fit
i, instantaneous Plot instantaneous growth rate vs time
h, help Show this message and exit.
The two tream instability simulation produces a field energy time trace that can be plotted with
pgkyl twostream_fieldEnergy.bp sel c0 pl logy
We have selected the \(E_x\) component because all other field components are essentially zero. We can then measure the growth rate in the field energy with the following command
pgkyl twostream_fieldEnergy.bp sel c0 growth
Such command produces the output below
fitGrowth: fitting region 100 > 7916
gamma = +3.79836e01 (current +5.814e02 R^2=7.407e01) 99.99% done [========= ]
gamma = +3.79836e01
This output indicates that it made the measurement taking into account all the data between the 100th and the 7916th point (if you examine the file twostram_fieldEnergy.bp you’d find that it has 7916 data points), and arrived at the growth rate \(\gamma=0.379836\), the the \(R^2\) of the fit being 0.7407.
The growth
command also allows for specifying a window in which to
perform the measurement via the minn
and maxn
flags. We could
then limit the window between the 2000th and the 6000th point with
pgkyl twostream_fieldEnergy.bp sel c0 growth minn 2000 maxn 6000
and the output would be
fitGrowth: fitting region 2000 > 6000
gamma = +3.79836e01 (current +4.322e01 R^2=9.998e01) 99.98% done [========= ]
gamma = +3.79836e01
giving the same result obtained above.
There is also an option for specifying a guess to \(A\) and \(\gamma\) in the fit, via the g flag:
pgkyl twostream_fieldEnergy.bp sel c0 growth minn 2000 maxn 6000 g 1. 0.36
fitGrowth: fitting region 2000 > 6000
gamma = +3.79836e01 (current +4.322e01 R^2=9.998e01) 99.98% done [========= ]
gamma = +3.79836e01
info¶
info
is a function of the postgkyl.data.Data
class that
returns information about the current dataset (see Data loading
for more details about the class itself). The function doesn’t take
any arguments and it is wrapped into the info
command.
Command line¶
Command help
pgkyl info help
Usage: pgkyl info [OPTIONS]
Print info of active datasets.
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
c, compact Show in compact mode
a, allsets All data sets.
h, help Show this message and exit.
info
always prints information about time stamp, number of
components (these are typically individual expansion coefficients
and/or vector components), number of dimensions with number of cells.
A useful capability of info
when working with complex pgkyl
commands (especially if using tags) is the ac
flags. For example,
if we had been working on the time averaged potential in the
integrate page, and instead of plotting we placed
info ac
at the end:
pgkyl gkionSound1x2vp1_phi_0.bp l '$t=0$' t phi0 \
"gkionSound1x2vp1_phi_[09]*.bp" t phis interp collect u phis t phiC \
integrate u phiC t phiInt 0 ev l 'Time average' t phiAvg 'phiInt 10. /' \
activate t phi0,phiAvg info ac
we would get the following output
Set $t=0$ (phi0#0)
Set 0 (phis#0)
Set 1 (phis#1)
...
Set 50 (phis#50)
Set collect (phiC#0)
Set (phiInt#0)
Set Time average (phiAvg#0)
Showing the label, tag and index of each dataset.
integrate¶
Integrate a dataset over specified direction(s).
Command line¶
Command help
$ pgkyl integrate help
Usage: pgkyl integr [OPTIONS] AXIS
Integrate data over a specified axis or axes
Options:
u, use TEXT Specify the tag to integrate
t, tag TEXT Optional tag for the resulting array
h, help Show this message and exit.
Consider the gyrokinetic simulation of an ion acoustic wave as an example. It outputs the integrated particle density over time, which can plot as follows:
pgkyl gkionSound1x2vp1_ion_intM0.bp pl
We can see from the values on the yaxis that the total number
of particles is 12.566. The number of particles should be conserved,
to machine precision. We can check this another way by integrating
the particle density along \(x\) (the 0th dimension) at the end
of the simulation with pgkyl
:
pgkyl gkionSound1x2vp1_ion_GkM0_1.bp interp integr 0 print
where we have abbreviated integrate
with integr
, and we use
the print command to print the result of the integral to screen. The
output of this command is simply
12.566370614359034
The integrate command can also be used to integrate higher dimensional datasets in one or more directions. We could take the ion distribution function and integrated along the \(v_\parallel\) and \(\mu\) directions (1st and 2nd dimensions, respectively) with
pgkyl gkionSound1x2vp1_ion_0.bp gkionSound1x2vp1_ion_GkM0_0.bp interp \
activate i 0 integr 1,2 ev l 'integrate 1,2' 'f[0] 6.283185 *' \
activate i1,2 pl f0 x 'x' y 'Number density, $n$'
Note
In this command we:
 First load the ion distribution function (*_ion_0.bp) and its number density (*_ion_GkM0_0.bp) at \(t=0\).
 Integrate the distribution function over velocity space with
activate i 0 integr 1,2
.  Multiply such integral by \(2\pi B_0/m_i\) (\(B_0=m_i=1\) here) with
ev l 'integrate 1,2' 'f[0] 6.283185 *'
.  Activate the number denstiy and integrated distribution function data sets and plot them with
activate i1,2 pl f0
.
and this should give approximately the same number density as the
GkM0
diagnostic outputted by the simulation, as shown below.
Another useful application of the integrate command is to integrate,
or average, over time (although note that the ev
command has a avg
operation that may make this easier). Usually
this requires collecting multiple frames into a single dataset with
the collect command, and then integrating over the
0th dimension (time).
So if we increase the tEnd
of the gyrokinetic ion sound wave
simulation to 10 and the number of frames to 50 we could
plot the electrostatic potential as a function of time and position
with
pgkyl "gkionSound1x2vp1_phi_[09]*.bp" interp collect pl x 'time' y 'x' clabel '$\phi$'
We can integrate this potential in time and plot it on top of the initial potential with
pgkyl gkionSound1x2vp1_phi_0.bp l '$t=0$' t phi0 \
"gkionSound1x2vp1_phi_[09]*.bp" t phis interp collect u phis t phiC \
integrate u phiC t phiInt 0 ev l 'Time average' t phiAvg 'phiInt 10. /' \
activate t phi0,phiAvg pl f0 x '$x$' y '$\phi$'
This command uses tags to select which dataset to perform an operation on. The end result is the plot below
showing that the time averaged potential is lower amplitude due to the collisionless Landau damping of the wave.
interpolate¶
Gkeyll’s simulation may employ higherorder methods with more than one degree of freedom per cell, typically expansion coefficients or nodal values. The interpolate command interpolates this higherorder representation onto a finer mesh with more points than the number of cells in the Gkeyll simulation.
Command line¶
Command help
pgkyl interpolate help
Usage: pgkyl interp [OPTIONS]
Interpolate DG data onto a uniform mesh.
Options:
b, basistype [msnsmomt] Specify DG basis.
p, polyorder INTEGER Specify polynomial order.
i, interp INTEGER Interpolation onto a general mesh of
specified amount.
u, use TEXT Specify a 'tag' to apply to (default all
tags).
t, tag TEXT Optional tag for the resulting array
r, read BOOLEAN Read from general interpolation file.
n, new for testing purposes
h, help Show this message and exit.
A Gkeyll simulation can have several degrees of freedom per cell. If we simply invoke the plot command we would then obtain a plot of all the coefficients in each cell. So for example, take the electron distribution function at the end of a two stream instability simulation and plot it with
pgkyl twostream_elc_100.bp pl x '$x$' y '$v_x$'
we would obtain figure below on the left, with all 8 DG expansion
coefficients plotted in phase space. Most commonly we are interested
in plotting the actual function that the expansion coefficients
represent, rather than plotting such coefficients. We can do that
by interpolating onto a finer mesh with the interpolate
command:
pgkyl twostream_elc_100.bp interp pl x '$x$' y '$v_x$'
which results in the figure below on the right.
We can compare the cell average of the electron distribution function with the interpolated distribution function with the following command
pgkyl twostream_elc_100.bp t fe interp u fe t fInterp sel u fe c0 t c0 \
ev l 'cell average' t cellAv 'c0 2 /' activ t fInterp,cellAv pl b
which divides the zeroth DG coefficient by 2 in order to obtain the cell average (1r 2D piecewise quadratic basis), and produces the following figure
Notice how the cell average (right) is naturally coarser grained, and the interpolated function (left) offers a smoother plot.
By default the interpolate
command interpolates onto a
uniform mesh that subdivides each cell in the simulation into
\(p+1\) cells in each direction, where \(p\) is the
polynomial order of the simulation. It is also possible to
interpolate onto finer meshes with the i
flag in order to
obtain smoother plots. However note that interpolating onto
finer meshes can also augment local maxima and/or minima. Below
we compare the final electron distribution function interpolated
onto a mesh with 3 subcells per cell in each direction (default
for \(p=2\)) and interpolating onto a mesh with 8 subcells
per cell in each direction:
pgkyl twostream_elc_100.bp t fe interp t i3 interp i 8 u fe t i8 activ t i3,i8 pl b x '$x$' y '$v_x$'
This example also shows the use of tags in order to tag datsets
and to instruct interpolate
which datasets to operate on
(via the u/use
flag). In order to request that
interpolate
operates on a given tagged dataset, one must pass
u
followed by the dataset we wish to interpolate. And in
order to create a new dataset outof the interpolated data
one must use the t
flag followed by the name (tag) of the
new dataset. In the above example the first interpolate
operates on the input data (no u
necessary because it
immediately precedes it and there is only one dataset at that point
in the chain) and creates a dataset tagged i3
. The second
interpolate
operates on the input data (u fe
) and creates
a dataset tagged i8
.
Script mode¶
interpolate
uses the GInterpModal
and GInterpNodal
classes based on the DG mode.
Parameter  Description  Default 

gdata (GData)  A GData object to be used.  
polyOrder (int)  The polynomial order of the discontinuous Galerkin discretization.  
basis (str)  The polynomial basis. Currently supported options are 'ns' for
nodal Serendipity, 'ms' for modal Serendipity, and 'mo'
for the maximal order basis. 
After the initialization, both GInterpModal
and GInterpNodal
can be used to interpolate data on a uniform grid and to calculate
derivatives
Member  Description 

interpolate(int component, bool stack) > narray, narray  Interpolates the selected component (default is 0) of the DG data on a uniform grid 
derivative(int component, bool stack) > narray, narray  Calculates the derivative of the DG data 
When the stack
parameter is set to true
(it is false
by
default), the grid and values are pushed to the GData
stack rather
than returned.
An example of the usage:
import postgkyl as pg
data = pg.data.GData('bgk_neut_0.bp')
interp = pg.data.GInterpModal(data, 2, 'ms')
iGrid, iValues = interp.interpolate()
listoutputs¶
You can list the unique filename stems of all the files written
by a simulation. This is useful if you are unsure about which
outputs were produced, or want to browse all possible files/diagnostics.
This command doesn’t take any arguments, and can be used with
the appreviation list
Command line¶
Command help
pgkyl list help
Usage: pgkyl list [OPTIONS]
List Gkeyll filename stems in the current directory
Options:
e, extension TEXT Output file extension [default: bp]
h, help Show this message and exit.
In the directory where the simulation was run type
pgkyl list
and you will obtain a list of the unique filename stems in that directory. For example, if I run the gkalfvenp1.lua simulation I would get
gkalfvenp1_allGeo
gkalfvenp1_apar
gkalfvenp1_dApardt
gkalfvenp1_electron
gkalfvenp1_electron_GkM0
gkalfvenp1_electron_GkM1
gkalfvenp1_electron_f0
gkalfvenp1_electron_f1
gkalfvenp1_ion
gkalfvenp1_ion_GkM0
gkalfvenp1_ion_GkM1
gkalfvenp1_laplacianWeight
gkalfvenp1_modifierWeight
gkalfvenp1_phi
Note that if multiple simulations with different input files have been run in the same directory it will list the unique filename stems for all such simulations.
mask¶
Command line¶
Command help
pgkyl mask h
Usage: pgkyl mask [OPTIONS]
Mask data with specified Gkeyll mask file.
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
f, filename TEXT Specify the file with a mask
l, lower FLOAT Specify the lower theshold to be masked out.
u, upper FLOAT Specify the upper theshold to be masked out.
h, help Show this message and exit.
Script mode¶
norm¶
Command line¶
Command help
pgkyl norm h
Usage: pgkyl norm [OPTIONS]
Normalize data
Options:
shift / noshift Shift minimal value to zero (default: False).
usefirst Normalize to first value in field
h, help Show this message and exit.
plot¶
postgkyl.output.plot
plots 1D and 2D data using Matplotlib.
The Python function is wrapped into the plot
command.
Contents
Default plotting¶
plot
automatically regnizes the dimensions of data and creates
either 1D line plot or 2D pcolormesh
plot using the Postgkyl
style file (Inferno color map).
Here is an example of 2D particle distribution function from the twostream instability simulation.
import postgkyl as pg
data = pg.Data('twostream_elc_80.bp')
dg = pg.GInterpModal(data)
dg.interpolate(stack=True)
pg.output.plot(data)
pgkyl twostream_elc_80.bp interpolate plot
Note that in this case the data does not contain the values of the distribution function directly but rather the expansion components of the basis functions. Therefore, interpolate was added to the flow to show the distribution function itself.
1D plots are created in a similar manner. For example, here is the electron density correfponding to the figure above.
import postgkyl as pg
data = pg.Data('twostream_elc_M0_80.bp')
dg = pg.GInterpModal(data)
dg.interpolate(stack=True)
pg.output.plot(data)
pgkyl twostream_elc_M0_80.bp interpolate plot
Plotting data with multiple components¶
Gkeyll data can contain multiple components. Typically, these are basis function expansion coefficients but can also correspond to components of a vector array like electromagnetic field or momentum. By default, Postgkyl plots each component into a separate subplot.
This can be seen if we do not use the interpolation from the previous example and let Postgkyl plot the expansion coefficients.
import postgkyl as pg
data = pg.Data('twostream_elc_M0_80.bp')
pg.output.plot(data)
pgkyl twostream_elc_M0_80.bp plot
Postgkyl automatically adds labels with component indices to each
subplot. If there are some labels already (either custom or when
working with multiple data sets), the component indices are
appended. Postgkyl also automatically calculates the numbers of rows
and columns (it tries to make a square). This can be overridden with
nSubplotRow
or nSubplotCol
.
The default behavior of putting each component to an individual
subplot can be supressed with the squeeze
parameter. This is
useful, for example, for comparing magnitudes. Note that the
magnitues of the expansion coefficients are quite different so this is
not the best example of the functionality.
import postgkyl as pg
data = pg.Data('twostream_elc_M0_80.bp')
pg.output.plot(data, squeeze=True)
pgkyl twostream_elc_M0_80.bp plot squeeze
Plotting multiple datasets¶
Postgkyl in a terminal can easily load multiple files (see Data loading for more details). By default, each data set creates its own figure.
pgkyl twostream_elc_70.bp twostream_elc_80.bp interp plot
Postgkyl automatically parses the names of the files and creates labels from the unique part of each one. Note that the labels can specified manually during Data loading.
This behavior can be supressed by specifying the figure to plot in. When the same figure is specified, data sets are plotted on top of each other.
import postgkyl as pg
data1 = pg.Data('twostream_elc_M0_70.bp')
dg = pg.GInterpModal(data1)
dg.interpolate(stack=True)
data2 = pg.Data('twostream_elc_M0_80.bp')
dg = pg.GInterpModal(data2)
dg.interpolate(stack=True)
pg.output.plot(data1, figure=0)
pg.output.plot(data2, figure=0)
pgkyl twostream_elc_M0_70.bp twostream_elc_M0_80.bp interp plot f0
Finally, the data sets can be added into subplots.
pgkyl twostream_elc_70.bp twostream_elc_80.bp interp plot f0 subplots
The same behavior can be achieved in a script as well but it requires slightly more manual control.
import postgkyl as pg
data1 = pg.Data('twostream_elc_M0_70.bp')
dg = pg.GInterpModal(data1)
dg.interpolate(stack=True)
data2 = pg.Data('twostream_elc_M0_80.bp')
dg = pg.GInterpModal(data2)
dg.interpolate(stack=True)
pg.output.plot(data1, figure=0, numAxes=2)
pg.output.plot(data2, figure=0, numAxes=2, startAxes=1)
Plotting modes¶
Appart from the default line 1D plots and continuous 2D plots, Postgkyl offers some additional modes.
Countour¶
import postgkyl as pg
data = pg.Data('twostream_elc_80.bp')
dg = pg.GInterpModal(data)
dg.interpolate(stack=True)
pg.output.plot(data, contour=True)
pgkyl twostream_elc_80.bp interpolate plot contour
Diverging¶
Diverging mode is similar to the default plotting mode but the colormap is changed to a redwhiteblue and the range is set to the plusminus maximum absolute value. It is particulary useful for visualizing changes, both in time and around a mean value.
Here we use the ev command to visualize the change from the initial conditions.
Group¶
In the group mode (maybe not the best name :/), one direction (either 0 or 1) is retained and the other is split into individual lineouts which are then plot over each other. The lines are colorcoded with the inferno colormap, i.e., from black to yellow as the coordinate increases. This could provide an additional insight into variation along one coordinate axis.
In the example, the 2D distribution function is first limited in the
first coordinate, z0
(in this case corresponding to x
), from
1.5 to 2.0 using the select command (otherwise there
would be too many lines). Then the plot with group=True
is used.
import postgkyl as pg
data = pg.Data('twostream_elc_80.bp')
dg = pg.GInterpModal(data)
dg.interpolate(stack=True)
pg.data.select(data, z0='1.5:2.0', stack=True)
pg.output.plot(data, group=1)
pgkyl twostream_elc_80.bp interpolate select z0 1.5:2.0 plot group 1
Formating¶
While Postgkyl is not necesarily meant for the production level figures for publications, it includes a decent amount of formating options.
The majority of a look of each figure, e.g., grid line style and
thickness or colormap, is set in a stule file. Custom matplotlib style
files can be specified with style
keyword. The default Postgkyl
style is the following:
figure.facecolor : white
lines.linewidth : 2
font.size : 12
axes.labelsize : large
axes.titlesize : 14
image.interpolation : none
image.cmap : inferno
image.origin : lower
grid.linewidth : 0.5
grid.linestyle : :
Labels¶
Postgkyl allows to specify all the axis labels and the plot title.
import postgkyl as pg
data = pg.Data('twostream_elc_80.bp')
dg = pg.GInterpModal(data)
dg.interpolate(stack=True)
pg.output.plot(data, xlabel=r'$x$', ylabel=r'$v_x$',
title=r'$k=\frac{1}{2}$')
pgkyl twostream_elc_80.bp interpolate \
plot xlabel '$x$' ylabel '$v_x$' \
title '$k=\frac{1}{2}$'
Axes and values¶
Postgkyl supports the logaritmic axes using the keywords logx
and logy
. In the example, the electric field energy is plotted
using the logarithmic yaxis to show the region of the linear growth
of the two stream instability. Note that Gkeyll stores \(E_x^2\),
\(E_y^2\), \(E_z^2\), \(B_x^2\), \(B_y^2\), and
\(B_z^2\) into six components of the fieldEnergy.bp
file. Therefore, the select command is used to plot only
the \(E_x^2\), which is the only component growing in this case.
import postgkyl as pg
data = pg.Data('twostream_fieldEnergy.bp')
pg.data.select(data, comp=0, stack=True)
pg.output.plot(data, logy=True)
pgkyl twostream_fieldEnergy.bp select c0 plot logy
Note that Postgkyl also incudes a diagnostic growth command that allows to fit the data with an exponential to get the growth rate.
Storing¶
Plotting outputs can be save as a PNG
files using the save
parameter which uses the data set name(s) to put together the name of
the image. Alternativelly, saveas
can be used to specify the
custom file name (without the extension, all the files are saved as
.png
). DPI of the result can be controlled with the dpi
parameter.
pr¶
Print data to screen.
Command line¶
Command help
pgkyl pr h
Usage: pgkyl pr [OPTIONS]
Print the data
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
h, help Show this message and exit.
This command allows us to print a dataset to screen. Large data sets like timetraces or multidimensional data might fill up the screen, so this might not be a good way to explore those. However, since postgkyl does not plot 0D data (single points), it is useful to print those to screen.
For example, we could integrate the initial number density in the two stream instability simulation and print out the total number of particles to the screen with
pgkyl twostream_elc_M0_0.bp interp integrate 0 pr
and produce the following output:
25.132741228817842
recovery¶
Command line¶
Command help
pgkyl recovery h
Usage: pgkyl recov [OPTIONS]
Interpolate DG data on a uniform mesh
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
t, tag TEXT Optional tag for the resulting array
b, basistype [msnsmo] Specify DG basis
p, polyorder INTEGER Specify polynomial order
i, interp INTEGER Number of poins to evaluate on
r, periodic Flag for periodic boundary conditions
c, c1 Enforce continuous first derivatives
h, help Show this message and exit.
Script mode¶
runchain¶
Command line usage¶
$ pgkyl runchain help
Usage: pgkyl runchain [OPTIONS]
Run the saved command chain
Options:
f, filename TEXT Specify file with stored chain (default
'pgkylchain.dat')
help Show this message and exit.
select¶
Select a subset of one or multiple datasets.
This subset can can be created by
 Choosing a cell along a direction (
z#
with an integer).  Choosing a range of cells along a direction (
z#
with a slice).  Choosing a specific component (
c
), which may be an expansion coefficient, fluid moment or, if used along with interpolate, a vector component. Multiple components can be selected with a slice of floats or with commaseparated values.  Evaluating the function at a specific coordinate in one direction
(using interpolate and
select z#
with a float) or segment in one dimension (using interpolate andselect z#
with a slize of floats).
Note
Postgkyl retains the number of dimensions and component index. This means that, for example, given a 16x16 2D dataset with 8 components per cell fixing the second coordinate and selecting one component will produce a dataset with a shape (16, 1, 1). For plotting purposes postgkyl treats such data as 1D.
Command line¶
Command help
pgkyl select help
Usage: pgkyl select [OPTIONS]
Subselect data from the active dataset(s). This command allows, for
example, to choose a specific component of a multicomponent dataset,
select a index or coordinate range. Index ranges can also be specified
using python slice notation (start:end:stride).
Options:
z0 TEXT Indices for 0th coord (either int, float, or slice)
z1 TEXT Indices for 1st coord (either int, float, or slice)
z2 TEXT Indices for 2nd coord (either int, float, or slice)
z3 TEXT Indices for 3rd coord (either int, float, or slice)
z4 TEXT Indices for 4th coord (either int, float, or slice)
z5 TEXT Indices for 5th coord (either int, float, or slice)
c, comp TEXT Indices for components (either int, slice, or coma
separated)
u, use TEXT Specify a 'tag' to apply to (default all tags).
t, tag TEXT Optional tag for the resulting array
h, help Show this message and exit.
Perhaps the two most common uses of select
are to choose a vector
component or to evaluate a function at a specific coordinate. Consider
the data produced by a
two stream instability VlasovMaxwell simulation.
The file twostream_field_0.bp
contains 24 components per cell.
That is because it contains 8 vector components (3 components of the
electric and magnetic fields each and 2 correction potentials, see
http://ammarhakim.org/sj/maxwelleigensystem.html) and 3
degrees of freedom per component (for piecewise linear basis in 1D).
If we wish to only see \(E_x\) then we can select it and plot it
with
pgkyl twostream_field_0.bp interp select c0 pl
which yields the following figure
We could also select multiple components with a commaseparated list
pgkyl twostream_field_0.bp interp select c0,3 pl x '$x$'
which yields the following plot of \(E_x\) and \(B_x\)
Or we could select all the components of the magnetic field in all frames with
pgkyl "twostream_field_[09]*.bp" interp select c3:6 pl f0 x '$x$' nolegend nsubplotrow 1
to show that this is an electrostatic simulation.
As a demonstration of using select
to evaluate functions at a given
coordinate we could evaluate the above \(E_x\) at \(x=\pi\), but
that would yield a single number. We could instead load all the field
data files, evaluate the \(x\)component of the electric field
at \(x=\pi\) in each frame, collect all those points into a single
dataset, and plot them as a function of time. This would be accomplished
with the following command
pgkyl "twostream_field_[09]*.bp" interp sel c0 z0 3.14159 collect pl x 'time' y '$E_x(x=\pi,t)$'
We will sometimes use the abbreviation sel
instead of select
for convenience. The resulting figure (above) shows how the amplitude of
the electric field at this point drastically increases as the instability
develops. There is an oscillatory behavior that is lost in this
exponential growth; if one zoomed into the earlier times of this plot,
such oscillations would become visible.
One can, and must, select coordinates at which to evaluate higher (than 2) dimensional datasets. Take the initial and final distribution functions of the electrons as an example; we can evaluate them at \(x=0\) and plot their variation along \(v_x\) with
pgkyl twostream_elc_0.bp l '$t=0$' twostream_elc_100.bp l '$t=50$' interp \
sel z0 0. pl x '$v_x$' y '$f_e(x=0,v_x)$' f0
producing the following plot:
If we were interested in investigating the distribution function in a specific cell along \(x\), say the first (0th in 0index), we could use
pgkyl twostream_elc_100.bp sel z0 0 interp pl x '$x$' y '$v_x$'
Notice that this produces a 2D colorplot because it takes the expansion coefficients in the 0th cell along \(x\) and all cells along \(v_x\), interpolates them onto a finer mesh and plots them (so the \(x\) extent of this plot is a single cell of the simulation). If we instead wanted a 1D plot of the distribution function along \(v_x\), we could first interpolate onto a finer mesh and then evaluate it at the 0th cell of the finer mesh using
pgkyl twostream_elc_100.bp interp sel z0 0 pl x '$x$' y '$f_e(x=6.25046,v_x)$'
or interpolate it and evaluate it at the lower boundary of the domain along \(x\), which is located at \(x=2\pi\), with
pgkyl twostream_elc_100.bp interp sel z0 6.283185 pl x '$x$' y '$f_e(x=2\pi,v_x)$'
These two commands are evaluating the distribution function at slightly different \(x\) coordinates (\(\Delta x/(p+1)/2\) apart to be precise, where \(\Delta x\) is the cell length of the simulation, and \(p\) the polynomial order of the basis). We can discern the difference between the two by plotting them together using the following command:
pgkyl twostream_elc_100.bp t fe interp sel t zfl z0 6.283185 sel u fe t zint z0 0 \
pl u zfl,zint f0 x '$v_x$' y '$f_e$'
This commmand used tags to indicate which dataset to perform the interpolation on, and to name the interpolated datasets. The result is
Lastly, we show that the select
command can also be used to restrict
and interpolated dataset to a segement along one direction when the
z#
flag is used with a slice of floats. For example, if one wants
to plot the initial electron distribution function at \(x=0\) for
positive velocities only, then one could employ
pgkyl twostream_elc_0.bp interp sel z0 0. z1 0.: pl x '$x$' y '$f_e(x=0,v_x,t=0)$'
Script mode¶
Parameter  Description  Default 

data (GData)  Data to subselect.  
coord0 (int, float, or slice)  Index corresponding to the first coordinate for the partial load. Either integer, float, or Python slice (e.g., ‘2:5’).  None 
coord1 (int, float, or slice)  Index corresponding to the second coordinate for the partial load. Either integer, float, or Python slice (e.g., ‘2:5’).  None 
coord2 (int, float, or slice)  Index corresponding to the third coordinate for the partial load. Either integer, float, or Python slice (e.g., ‘2:5’).  None 
coord3 (int, float, or slice)  Index corresponding to the fourth coordinate for the partial load. Either integer, float, or Python slice (e.g., ‘2:5’).  None 
coord4 (int, float, or slice)  Index corresponding to the fifth coordinate for the partial load. Either integer, float, or Python slice (e.g., ‘2:5’).  None 
coord5 (int, float, or slice)  Index corresponding to the sixth coordinate for the partial load. Either integer, float, or Python slice (e.g., ‘2:5’).  None 
comp (int, slice, or multiple)  Index corresponding to the component for the partial load. Either integer, Python slice (e.g., ‘2:5’), or multiple.  None 
Unlike for the partial load parameters (see :ref:), float point numbers can be specified instead of just integers. In that case, Postgkyl treats it as a grid value and automatically finds and index of a grid point with the closest value. This works both for the single index and for specifying a slice.
In order to use select in a Python script one must interpolate the
nodal/modal dataset inplace (normally it produces a new dataset, i.e. outofplace)
and pass the original Gkeyll data to the select
command. For example, in order to
select the 0th coordinate at the value 0.0 we would use:
import postgkyl as pg
pgData = pg.GData(fileName)
pgInterp = pg.GInterpModal(pgData, polyOrder, basisType)
pgInterp.interpolate(overwrite=True)
x, f_z0eq0p0 = pg.data.select(pgData, z0=0.0)
where pg.data.select
returns the new grid (x) and the new field (f_z0eq0p0).
tenmoment¶
Command line¶
Command help
pgkyl tenmoment help
Usage: pgkyl tenmom [OPTIONS]
Extract tenmoment primitive variables from tenmoment conserved
variables.
Options:
u, use TEXT Specify a 'tag' to apply to (default all
tags).
v, variable_name [densityxvelyvelzvelvelpressureTensorpxxpxypxzpyypyzpzzpressure]
Variable to plot
h, help Show this message and exit.
Script mode¶
trajectory¶
Command line¶
Command help
pgkyl trajectory h
Usage: pgkyl traj [OPTIONS]
Animate a particle trajectory
Options:
fixaspect Enforce the same scaling on both axes.
show / noshow Turn showing of the plot ON and OFF (default:
ON).
i, interval INTEGER Specify the animation interval.
save Save figure as PNG.
velocity / novelocity Plot velocity vectors.
saveas TEXT Name to save the plot as.
e, elevation FLOAT Set elevation
a, azimuth FLOAT Set azimuth
n, numframes INTEGER Set number of frames for the animation
xmin FLOAT Minimum value of the xcoordinate
xmax FLOAT Maximum value of the xcoordinate
ymin FLOAT Minimum value of the ycoordinate
ymax FLOAT Maximum value of the ycoordinate
zmin FLOAT Minimum value of the zcoordinate
zmax FLOAT Maximum value of the zcoordinate
u, use TEXT Specify a 'tag' to apply to (default all tags).
h, help Show this message and exit.
val2coord¶
Command line¶
Command help
pgkyl val2coord h
Usage: pgkyl val2 [OPTIONS]
Given a dataset (typically a DynVector) selects columns from it to create
new datasets. For example, you can choose say column 1 to be the Xaxis of
the new dataset and column 2 to be the Yaxis. Multiple columns can be
choosen using range specifiers and as many datasets are then created.
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
t, tag TEXT Tag for the result
x TEXT Select components that will became the grid of the new
dataset.
y TEXT Select components that will became the values of the new
dataset.
h, help Show this message and exit.
velocity¶
Command line¶
Command help
pgkyl velocity h
Usage: pgkyl vel [OPTIONS]
Options:
d, density TEXT Tag for density [default: density]
m, momentum TEXT Tag for momentum [default: momentum]
t, tag TEXT Tag for the result [default: velocity]
l, label TEXT Custom label for the result [default: velocity]
h, help Show this message and exit.
write¶
Postgkyl can store data into a new ADIOS bp
file (default; useful
for storing partially processed data) or into a ASCII txt
file
(useful when one wants to open the data in a program that does not
support bp
or h5
).
Command line usage¶
Command help
pgkyl write help
Usage: pgkyl writ [OPTIONS]
Write active dataset to a file. The output file format can be set with
``mode``, and is ADIOS BP by default. If data is saved as BP file it can
be later loaded back into pgkyl to further manipulate or plot it.
Options:
u, use TEXT Specify a 'tag' to apply to (default all tags).
f, filename TEXT Output file name
m, mode [bptxtnpy] Output file mode. One of `bp` (ADIOS BP file;
default), `txt` (ASCII text file), or `npy` (NumPy
binary file)
b, buffersize INTEGER Set the buffer size for ADIOS write (default: 1000
MB)
h, help Show this message and exit.
Script mode¶
The write
command internally calls the write()
method of the
GData
class.
import postgkyl as pg
data = pg.data.GData('bgk_neut_0.bp')
data.write()
Publications and theses¶
Note
We acknowledge funding from the Department of Energy SciDAC programs (MGK, PI David Hatch and HPBS, PI CS Chang) as well as National Science Foundation via various grants to PPPL and Virgina Tech, the AirForce Office of Scientific Research, National Aeronautics and Space Administration, DOE’s APRAE program, graduate student fellowships via DOE’s CSGF and NASA’s FINESST programs, and postdoctoral fellowships via DOE’s ORISE and NSF AGS PRF programs.
A good source of various benchmarks and other tests is A. Hakim’s Simulation Journal and its github webpage.
We have also compiled input files for the simulations reported in publications in this repository. Note that this collection is incomplete as not all authors have desposited they input files with us.
Doctoral Dissertations¶
 Mandell, N. R. (2021, March 26) “Magnetic Fluctuations in Gyrokinetic Simulations of ScrapeOff Layer Turbulence”. Ph.D. dissertation, Princeton University, 2021. arXiv:2103.16062
 Juno, J. (2020, March 27) “A Deep Dive into the Distribution Function: Understanding Phase Space Dynamics Using Continuum VlasovMaxwell Simulations”. Ph. D. dissertation, University of Maryland, College Park, 2020. arXiv:2005.13539
 Bernard, T. N. “Discontinuous Galerkin Modeling of Plasma Turbulence in a Simple Magnetized Torus”. Ph. D. dissertation, The University of Texas at Austin, 2019. PDF
 Ng, J. “Fluid closures for the modeling of reconnection and instabilities in magnetotail current sheets”. Ph.D. dissertation, Princeton University, 2018. PDF
 Cagas, P. (2018, July 30). “Continuum kinetic simulations of plasma sheaths and instabilities”. Ph.D. dissertation, Virginia Polytechnic Institute and State University, 2018. https://vtechworks.lib.vt.edu/handle/10919/84979
 Shi, E. L. (2017, August 24). “Gyrokinetic Continuum Simulation of Turbulence in OpenFieldLine Plasmas”, Ph.D. dissertation, Princeton University, 2017 arXiv:1708.07283
Algorithms papers¶
 Francisquez, M., Mandell, N. R., Hakim, A., Hammett, G. W. (2021) “Mapped discontinuous Galerkin interpolations and sheared boundary conditions” arXiv2110.02249
 Cagas, P and Hakim, A and Srinivasan, B. (2021) “A boundary value “reservoir problem” and boundary conditions for multimoment multifluid simulations of sheaths.” Physics of Plasmas 28.1. https://doi.org/10.1063/5.0024510
 Hakim, A and Juno, J. (2020). “Aliasfree, Matrixfree, and Quadraturefree Discontinuous Galerkin Algorithms for (Plasma) Kinetic Equations”. SC20: Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis, IEEE Press. https://doi.org/10.1109/SC41405.2020.00077
 Francisquez, M., Bernard, T. N., Mandell, N. R., Hammett, G. W., Hakim, A. (2020). “Conservative discontinuous Galerkin scheme of a gyroaveraged Dougherty collision operator”, Nuclear Fusion, 60, (9). https://doi.org/10.1088%2F17414326%2Faba0c9
 Hakim, A., Francisquez, M., Juno, J., & Hammett, G. W. (2020). “Conservative discontinuous Galerkin schemes for nonlinear Dougherty–Fokker–Planck collision operators”, Journal of Plasma Physics, 86, (4). https://doi.org/10.1017/S0022377820000586
 Wang, L., Hakim, A., Ng, J., Dong, C., & Germaschewski, K. (2020). “Exact and locally implicit source term solvers for multifluidMaxwell systems”, Journal of Computational Physics, 415, 109510. https://doi.org/10.1016/j.jcp.2020.109510
 Cagas, P., Hakim, A., & Srinivasan, B. (2020). “Plasmamaterial boundary conditions for discontinuous Galerkin continuumkinetic simulations, with a focus on secondary electron emission”, Journal of Computational Physics, 406, 109215. https://doi.org/10.1016/j.jcp.2019.109215
 Mandell, N. R., Hakim, A., Hammett, G. W., & Francisquez, M. (2020). “Electromagnetic fullf gyrokinetics in the tokamak edge with discontinuous Galerkin methods”, Journal of Plasma Physics, 86. https://doi.org/10.1017/S0022377820000070
 Juno, J., Hakim, A., TenBarge, J., Shi, E., & Dorland, W. (2018). “Discontinuous Galerkin algorithms for fully kinetic plasmas”, Journal of Computational Physics, 353, 110–147. https://doi.org/10.1016/j.jcp.2017.10.009
 Hakim, A., Hammett, G. W., Shi, E. L. (2014). “On discontinuous Galerkin discretizations of secondorder derivatives”, arXiv:1405.5907
Physics papers¶
 Mandell, N. R., Hammett, G. W., Hakim, A., Francisquez, M. (2021). “Reduction of transport due to magnetic shear in gyrokinetic simulations of the scrapeoff layer”, arXiv:2112.14220.
 Francisquez, M., Juno, J., Hakim, A., Hammett, G. W., Ernst, D. R. (2021). “Improved multispecies Dougherty collisions”, arXiv:2109.10381.
 Mandell, N. R., Hammett, G. W., Hakim, A., Francisquez, M. (2021). “Turbulent broadening of electron heatflux width in electromagnetic gyrokinetic simulations of a helical scrapeoff layer model”, arXiv:2112.06880.
 Mathews, A., Mandell, M., Francisquez, M., Hughes, J., Hakim, A. (2021). “Turbulent field fluctuations in gyrokinetic and fluid plasmas”, arXiv:2107.09744. Accepted for publication in Physics of Plasmas.
 Wang, L., Hakim, A., Srinivasan, B., Juno, J. (2021). “Electron cyclotron drift instability and anomalous transport: twofluid moment theory and modeling”, arXiv:2107.09874.
 Jenab, S. M., Brodin, G., Juno, J., Kourakis, I. (2021). “Ultrafast Electron Holes in Plasma Phase Space Dynamics”, Scientific Reports, 11. https://doi.org/10.1038/s4159802195652w
 Juno, J., Howes, G. G., TenBarge, J. M., Wilson III, L. B., Spitkovsky, A., Caprioli, D., Klein, K. G., Hakim, A. (2021). “A fieldparticle correlation analysis of a perpendicular magnetized collisionless shock”, Journal of Plasma Physics, 87, (3). https://doi.org/10.1017/S0022377821000623
 Pezzi, O., Liang, H., Juno, J. L., Cassak, P. A., Vásconez, C. L., SorrisoValvo, L., Perrone, D., Servidio, S., Roytershteyn, V., TenBarge, J.M., & Matthaeus, W. H. (2021). “Dissipation measures in weakly collisional plasmas”, Monthly Notices of the Royal Astronomical Society, 505, (4), Pages 4857–4873. https://doi.org/10.1093/mnras/stab1516
 Ng, J., Hakim, A., Wang, L., & Bhattacharjee, A. (2020). “An improved tenmoment closure for reconnection and instabilities”, Physics of Plasmas, 27, 082106. https://doi.org/10.1063/5.0012067
 Juno, J., Swisdak, M. M., TenBarge. J. M., Skoutnev, V., & Hakim, A. (2020). “Noiseinduced magnetic field saturation in kinetic simulations”, Journal of Plasma Physics, 86, (4). https://doi.org/10.1017/S0022377820000707
 Ng, J., Chen, L.‐J., Hakim, A., & Bhattacharjee, A. (2020). “Reconstruction of electron and ion distribution functions in a magnetotail reconnection diffusion region”, Journal of Geophysical Research: Space Physics, 125, e2020JA027879. https://doi.org/10.1029/2020JA027879
 Francisquez, M., Bernard, T. N., Zhu, B., Hakim, A., Rogers, B. N., & Hammett, G. W. (2020). “Fluid and gyrokinetic turbulence in open fieldline, helical plasmas”, Physics of Plasmas, 27, 082301. https://doi.org/10.1063/5.0005333
 Bernard, T. N., StoltzfusDueck, T., Gentle, K. W., Hakim, A., Hammett, G. W., & Shi, E. L. (2020). “Investigating shear flow through continuum gyrokinetic simulations of limiter biasing in the Texas Helimak”, Physics of Plasmas, 27, 062304. https://doi.org/10.1063/5.0003904
 Hakim, A. H., Mandell, N. R., Bernard, T. N., Francisquez, M., Hammett, G. W., & Shi, E. L. “Continuum electromagnetic gyrokinetic simulations of turbulence in the tokamak scrapeoff layer and laboratory devices”, Physics of Plasmas, 27, 042304. https://doi.org/10.1063/1.5141157
 Pusztai, I., Juno, J., Brandenburg, A., Tenbarge, J. M., Hakim, A., Francisquez, M., & Sundström, A. (2020). “Dynamo in Weakly Collisional Nonmagnetized Plasmas Impeded by Landau Damping of Magnetic Fields”, Physical Review Letters, 124, 255102. https://link.aps.org/doi/10.1103/PhysRevLett.124.255102
 TenBarge, J. M., Ng, J., Juno, J., Wang, L., Hakim, A. & Bhattacharjee, A. (2019). “An extended MHD study of the 16 October 2015 MMS diffusion region crossing”, Journal of Geophysical Research: Space Physics, 124, 84748487. https://doi.org/10.1029/2019JA026731
 Ng, J., Hakim, A., Juno, J., & Bhattacharjee, A. (2019). Drift instabilities in thin current sheets using a two‐fluid model with pressure tensor effects. Journal of Geophysical Research: Space Physics, 124, 33313346. https://doi.org/10.1029/2018JA026313
 Dong, C., Wang, L., Hakim, A., Bhattacharjee, A., Slavin, J. A., DiBraccio, G. A., & Germaschewski, K. (2019). “A Novel TenMoment Multifluid Model for Mercury: From the Planetary Conducting Core to the Dynamic Magnetosphere”, Geophysical Review Letters, 46, 1158411596. https://doi.org/10.1029/2019GL083180
 Shi, E. L., Hammett, G. W., StoltzfusDueck, T., & Hakim, A. (2019). “Fullf gyrokinetic simulation of turbulence in a helical openfieldline plasma”, Physics of Plasmas, 26, 012307. https://doi.org/10.1063/1.5074179
 Bernard, T. N., Shi, E. L., Gentle, K. W., Hakim, A., Hammett, G. W., StoltzfusDueck, T., & Taylor, E. I. (2019). “Gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak”, Physics of Plasmas, 26, 042301. https://doi.org/10.1063/1.5085457
 Skoutnev, V., Hakim, A., Juno, J., & TenBarge, J. M. (2019). “TemperatureDependent Saturation of WeibelType Instabilities in Counterstreaming Plasmas”, Astrophysical Journal Letters, 872, (2). https://doi.org/10.3847%2F20418213%2Fab0556
 Sundström, A., Juno, J., TenBarge, J. M., & Pusztai, I. (2019). “Effect of a weak ion collisionality on the dynamics of kinetic electrostatic shocks”, Journal of Plasma Physics, 85. https://doi.org/10.1017/S0022377819000023
 Srinivasan, B. and Hakim, A. (2018). “Role of electron inertia and electron/ion finite Larmor radius effects in lowbeta, magnetoRayleighTaylor instability”, Physics of Plasmas, 25, 092108. https://doi.org/10.1063/1.5046098
 Ng, J., Hakim, A., & Bhattacharjee, A. (2018). “Using the maximum entropy distribution to describe electrons in reconnecting current sheets”, Physics of Plasmas, 25, 082113. https://doi.org/10.1063/1.5041758
 Wang, L., Germaschewski, K., Hakim, A., Dong, C., Raeder, J., & Bhattacharjee, A. (2018). “Electron Physics in 3D TwoFluid 10Moment Modeling of Ganymede’s Magnetosphere”, Journal of Geophysical Research: Space Physics, 41 (A3), 8688–16. https://doi.org/10.1002/2017JA024761
 Pusztai, I., TenBarge, J. M., Csapó, A. N., Juno, J., Hakim, A., Yi, K & Fülöp, T. (2018). “Low Machnumber collisionless electrostatic shocks and associated ion acceleration”, Plasma Physics and Controlled Fusion, 60 (3), 035004–11. https://doi.org/10.1088/13616587/aaa2cc
 Shi, E. L., Hammett, G. W., StolzfusDueck, T., Hakim, A. (2017). “Gyrokinetic continuum simulation of turbulence in a straight openfieldline plasma”, Journal of Plasma Physics, 83, 1–27. https://doi.org/10.1017/S002237781700037X
 Cagas, P., Hakim, A., Scales, W., Srinivasan, B. (2017). “Nonlinear saturation of the Weibel instability”, Physics of Plasmas, 24 (11), 112116. https://doi.org/10.1063/1.4994682
 Ng, J., Hakim, A., Bhattacharjee, A., Stanier, A., & Daughton, W. (2017). “Simulations of antiparallel reconnection using a nonlocal heat flux closure”, Physics of Plasmas, 24 (8), 082112. https://doi.org/10.1063/1.4993195
 Stanier, A., Daughton, W., Simakov, A. N., Chacón, L., Le, A., Karimabadi, H., Ng, J., & Bhattacharjee, A. (2017). “The role of guide field in magnetic reconnection driven by island coalescence”, Physics of Plasmas, 24, 022124. https://doi.org/10.1063/1.4976712
 Cagas, P., Hakim, A., Juno, J., Srinivasan, B. (2017). “Continuum kinetic and multifluid simulations of classical sheaths”, Physics of Plasmas, 24 (2), 022118. https://doi.org/10.1063/1.4976544
 Ng, J., Huang, Y.M., Hakim, A., Bhattacharjee, A., Stanier, A., Daughton, W., Wang, L., & Germaschewski, K. (2015). “The island coalescence problem: Scaling of reconnection in extended fluid models including higherorder moments”, Physics of Plasma, 22, 112104. https://doi.org/10.1063/1.4935302
 Stanier, A., Daughton, W., Chacón, L., Karimabadi, H., Ng, J., Huang, Y.M., Hakim, A., & Bhattacharjee, A. (2015). “Role of Ion Kinetic Physics in the Interaction of Magnetic Flux Ropes”, Physical Review Letters, 115, 175004. https://doi.org/10.1103/PhysRevLett.115.175004
 Wang, L., Hakim, A. H., Bhattacharjee, A., & Germaschewski, K. (2015). “Comparison of multifluid moment models with particleincell simulations of collisionless magnetic reconnection”, Physics of Plasmas, 22 (1), 012108. https://doi.org/10.1063/1.4906063
Presentations¶
You can browse a folder of pdf / PowerPoint / Keynote files of Gkeyll presentations, or click on links below.
Note
Google Drive often does not display the slides properly. Please download the slides on your local machine to view them properly.
2021
 “Prioritizing, Leveraging, and Disseminating Fundamental Algorithms Research”, Jimmy Juno, PPPL CSD Seminar, November 2021. pdf
 “Modelling AUG scrapeofflayer plasma with fullf continuum Electromagnetic Gyrokinetic simulation”, Rupak Mukherjee, APS DPP Annual Meeting, November 2021. pdf, Keynote (with movies), video recording of talk
 “Electromagnetic fullf continuum gyrokinetic simulation of plasma turbulence in scrapeoff layer of ASDEX Upgrade”, Rupak Mukherjee, 19th European Fusion Theory Conference (Virtual in Consorzio RFX), October 2021. pdf, Keynote (with movies)
 “Electromagnetic fullf continuum gyrokinetic simulation of plasma turbulence in scrapeoff layer of ASDEX Upgrade tokamak”, Rupak Mukherjee, 5th Asia Pacific Conference on Plasma Physics (Virtual), September 2021. pdf, Keynote (with movies), video recording of talk
 “Modelling AUG scrapeofflayer plasma with fullf continuum Electromagnetic Gyrokinetic simulation”, Rupak Mukherjee, Virtual 25th Joint EUUS TTF Meeting, September 2021. pdf
 “Simulation of AUG scrapeofflayer plasma with fullf continuum Electromagnetic Gyrokinetic simulation”, Rupak Mukherjee, Sherwood Fusion Theory conference, August 2021. pdf, Keynote (with movies), video recording of talk
 “Initial Gkeyll simulations of ScrapeOffLayer Turbulence in ASDEXU”, Rupak Mukherjee, MPPC Meeting, January 2021. pdf, Keynote (with movies), video recording of talk
2020
 “Electromagnetic fullf gyrokinetic simulation of ASDEX SOL turbulence with discontinuous Galerkin method”, Rupak Mukherjee, APS DPP Annual Meeting, November 2020. pdf, Keynote (with movies), video recording of talk
 “Investigating magnetic fluctuations in gyrokinetic simulations of tokamak SOL turbulence”, Noah Mandell, APS DPP Annual Meeting, November 2020. pdf, Keynote (with movies), video recording of talk
 “A Deep Dive into the Distribution Function: Understanding Phase Space Dynamics Using Continuum VlasovMaxwell Simulations”, Jimmy Juno, APS DPP Annual Meeting, November 2020. Keynote
 “Balancing Flexibility and Usability in the Gkeyll Simulation Framework”, Jimmy Juno, APS DPP Annual Meeting, November 2020 Keynote
 “Studies of plasma sheaths using novel numerical schemes with selfconsistent emitting walls and FokkerPlanck collisions”, Petr Cagas, APS DPP Annual Meeting, November 2020. pdf, video recording
 “Kinetic Boltzmann modeling of neutral transport for a continuum gyrokinetic code”, Tess Bernard, APS DPP Annual Meeting, November 2020. pdf
 “Aliasfree, Matrixfree, and Quadraturefree Discontinuous Galerkin Algorithms for (Plasma) Kinetic Equations”, Ammar Hakim. SuperComputing 2020, November 2020. ppt
 “Investigating magnetic fluctuations in gyrokinetic simulations of tokamak SOL turbulence”, Noah Mandell, PPPL Theory Research & Review Seminar, October 2020. pdf, Keynote (with movies)
 “Investigating magnetic fluctuations in tokamak SOL turbulence using Gkeyll gyrokinetic simulations”, Noah Mandell, PPPL Monthly Research Meeting, October 2020. pdf, Keynote (with movies)
 “Magnetic fluctuations in gyrokinetic simulations of tokamak SOL turbulence”, Noah Mandell, Journal of Plasma Physics Frontiers colloquium series, April 2020. pdf, Keynote (with movies)
 “Initial SOL turbulence results from the Gkeyll code, including first electromagnetic effects”, Greg Hammett, AUG Seminar, Garching, January 2020. pdf, ppt (with movies)
2019
 “Continuum Electromagnetic Gyrokinetic Simulations of Turbulence in the Tokamak ScrapeOff Layer and Laboratory Devices”, Ammar Hakim, APS Division of Plasma Physics, Fort Lauderdale, 2019.
 “Tests of a Discontinuous Galerkin scheme for Hamiltonian systems in noncanonical coordinates”, Rupak Mukherjee, APS Division of Plasma Physics, Fort Lauderdale, 2019. pdf
 “Gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak”, Tess Bernard, Sherwood Fusion Theory Conference, Princeton, April 2019.
 “Gyrokinetic continuum simulations of plasma turbulence in the Texas Helimak”, Tess Bernard, 24th Joint USEU Transport Task Force Meeting, Austin, March 2019.
2016
 “FullF gyrokinetic simulations of the LAPD device with open field lines and sheath boundary conditions”, Greg W. Hammett, Eric L. Shi, Ammar Hakim, Oxford Plasma Theory Group Seminar, Nov. 17, 2016. pdf, ppt
Not very complete. more to be added…
Developer notes¶
On use of the Maxima CAS¶
Throughout Gkyl a large amount of Lua and C++ code is automatically pregenerated using the Maxima computer algebra system (CAS). Maxima is free and has a vast amount of features. For some of the calculations the use of a CAS is essential as the algebra, even though relatively easy, is very tedious, needing thousands of evaluations of various integrals etc.
A very pleasant frontend for Maxima is provided by the wxmaxima program. This provides a “document based” interface to Maxima and one can mix regular text (and equations) with Maxima interactions.
A very comprehensive physics oriented tutorial is Maxima by Example by Edwin Woollett.
All Maxima code is checked into the gkyl/casscripts directory. To use the Maxima code in this directory you need to tell Maxima to find it. To do this, create the directory (if it does not exist already):
mkdir $HOME/.maxima
In this create or edit the file called “maximainit.mac” and add the following lines to it:
file_search_maxima: append(file_search_maxima,
["PATH_TO_YOUR_GKYL/gkyl/casscripts/###.{lisp,mac}"]) $
Where “PATH_TO_YOUR_GKYL” is the full path to the location where your
gkyl source lives. Start/restart Maxima. Once you do this, then the
Maxima code in the casscripts
directory can be loaded, for
example as:
load("modalbasis")$
load("basisprecalc/basisSer1x1v")$
This will load the code to work with Modal basis functions and the serendipity basis sets in 1x1v into your Maxima session/code.
To make plots on Maxima, you can use the excellent draw2d/3d packages. Chapter 4 of this manual describes the draw packages. To get plotting to work you need to install Gnuplot and set some paths properly. On a Mac, the maximainit.mac file looks like:
load("draw")$
gnuplot_command: "/Applications/Gnuplot.app/Contents/Resources/bin/gnuplot" $
set_plot_option([gnuplot_term, "qt"],
[gnuplot_preamble, "set object rectangle from screen 0,0 to screen 1,1 behind fillcolor rgb 'white' fillstyle solid noborder"]
)$
set_draw_defaults(terminal=qt,
user_preamble="set object rectangle from screen 0,0 to screen 1,1 behind fillcolor rgb 'white' fillstyle solid noborder",
nticks=200,
line_width=2
)$
file_search_maxima: append(file_search_maxima,
["PATH_TO_YOUR_GKYL/gkyl/casscripts/###.{lisp,mac}"]) $
Again, remember “PATH_TO_YOUR_GKYL” is the full path to the location where your gkyl source lives. On Linux or Windows you will need to experiment with paths and settings to get plots to work.
Note to developers on maximum default available memory with sbcl
LISP compiler¶
Depending on which LISP compiler Maxima is using, by default, the compiler may not be able to claim/use all of the RAM on your computer. To change this a developer needs to edit the “maxima” executable themselves. On MacOS go to the following directory (assuming you have installed Maxima in your Applications):
cd /Applications/Maxima.app/Contents/Resources/opt/bin
Once in this directory, with your favorite text editor, open up the file maxima
.
Search for the specific LISP compiler being used. For example, the most common LISP compiler
amongst the Gkeyll development team is sbcl
. You should see a conditional statement like:
elif [ "$MAXIMA_LISP" = "sbcl" ]; then
within the conditional add:
MAXIMA_LISP_OPTIONS+="dynamicspacesize 100000"
or a larger number depending on the size of your RAM.
Modal basis functions¶
Gkyl uses orthonormal modal basis functions. The basis functions are defined on a ddimensional hypercube \(I_d = [1,1]^d\). Let \(\psi_k(\mathbf{x})\), \(k=1,\ldots,N\), with \(\psi_1\) a constant, be the basis. Then the basis satisfy
where the angle brackets indicate integration over the hypercube \(I_d\). In this note I describe some common operations that are needed while working with these basis sets. (All of this relatively straightforward stuff, but it is good to write it down somewhere).
Contents
Precomputed basis functions¶
Precomputed modal basis sets on phasespace for 1x1v, 1x2v, 1x3v, 2x2v, 2x3v, 3x2v, and 3x3v phasespace are stored as Lisp files in gkyl/casscripts/basisprecalc directory. Both maximalorder and serendipity basis sets are computed for polynomial orders 1, 2, 3 and 4. Computing orthonormal basis set in higher dimensions is timeconsuming and so these precomputed lisp files should be used. For example:
load("basisprecalc/basisSer2x3v")
will load the serendipity basis sets in 2x3v phasespace.
NOTE: please read Maxima notes to get this command to work.
The files define the following Maxima variables:
varsC, varsP, basisC, basisP
The varsC are the configurationspace independent variables (\(x,y\) in 2X) and varsP are the phasespace independent variables (\(x,y,vx,vy,vz\) in 2X3V). The ith entry in the list basisC and basisP contain the basis sets of polynomial order \(i\). Hence, to access the phasespace basis for \(p=2\) do something like:
basisP2 : basisP[2].
Computing cell averages¶
Consider some function that is expanded in the basis:
where \(\boldsymbol{\eta}(\mathbf{x})\) maps the physical space to logical space. Then, the cellaverage is defined as
where \(\psi_1\) is constant. By orthonormality we have \(\langle \psi_1^2 \rangle = 1\) which indicates that
This means that the cellaverage is given by
Convolution of two functions¶
Now consider two functions, \(f\) and \(g\) that are both expanded in the basis. The inner product of the two functions is:
Note the the bar on the LHS, i.e. this expression gives the average of the inner product on a single cell. Now, as the basis functions are orthonormal, this leads to the particularly simple expression
This makes it very easy to compute things like electromagnetic energy, and other quadratic quantities.
The recovery Maxima code¶
A major algorithm used in many parts of Gkeyll is the recovery procedure. This is used to construct DG (hyper)diffusion equation solvers, the LBO diffusion kernels and even some very highorder hyperbolic solvers. We have attempted to abstract the recovery procedure (which is rather complicated) into three very simple to use functions.
StrongStability preserving RungeKutta timesteppers¶
The Gkyl DG solvers use SSPRK timesteppers. Three steppers are implemented: SSPRK2, SSPRK3 and a fourstage SSPRK3 that allows twice the CFL (for the cost of additional memory) as the other schemes. See [DurranBook] page 56. The schemes are described below. Here, the symbol \(\mathcal{F}\) is used to indicate a firstorder Euler update:
where \(\mathcal{L}[f]\) is the RHS operator from the spatial discretization of the DG scheme.
Contents
SSPRK2¶
with \(CFL \le 1\).
SSPRK3¶
with \(CFL \le 1\). As this scheme has three stages instead of two, it will take about \(1.5X\) longer to run than the SSPRK2 scheme.
Four stage SSPRK3¶
with \(CFL\le 2\). Note that this scheme has four stages, but allows twice the timestep that SSPRK2 and SSPRK3, hence will result in a speed up of \(1.5X\) compared to the threestage SSPRK3 scheme.
Region of absolute stability¶
For each of the above schemes, I have plotted below the region of absolute stability. Note that only the RK3 schemes are stable when there is no diffusion in the system, and hence should be prefered.
References¶
[DurranBook]  Dale E. Durran, “Numerical Methods for Fluid Dynamics”, Springer. Second Edition. 
Normalized units for the VlasovMaxwell system¶
The Gkyl design philosophy involves the implementation of unitfull systems of equations, i.e., Gkyl simulations can be run with real parameters for direct comparison with experiments, with universal constants specified by using values provided by the National Institute of Standards and Technology. For example, in the absence of collisions, the VlasovMaxwell system of equations in S.I. units is as follows,
However, one may not always wish to run simulations with the unitfull system. Instead, one can consider a normalized set of equations. A natural choice for the normalization of the VlasovMaxwell system of equations would redefine all the relevant quantities as follows,
where
are the speed of light, electron plasma frequency, and the electron skin depth respectively. Note that the charge normalization means that in a protonelectron plasma, \(\tilde{q}_s = \pm 1\), and the density normalization is such that in a quasineutral plasma, the initial density of each species is 1.0. We can also check that the electric and magnetic field normalizations are reasonable by making sure that the normalization has the correct units for the electric and magnetic fields in S.I. units,
With these normalizations, the VlasovMaxwell system of equations then becomes,
This system of equations has the obvious advantage that universal constants, such as \(\epsilon_0\), are eliminated. In doing so, one does not need to worry about the propagation of round off error from, for example, the accumulation of the current to the electric field in the AmpereMaxwell law, \(E^{n+1} = E^{n} + \delta t \mathbf{J}/\epsilon_0\) becomed \(E^{n+1} = E^{n} + \delta t \tilde{\mathbf{J}}\). Given Gkyl’s unitfull representation, a simple way to force the VlasovMaxwell solver to “use” these units is to specify the following parameters be equal to 1.0,
With the above parameters set to 1.0, then VlasovMaxwell simulations require only a few parameters to be completely determined. In a protonelectron plasma, these are, the protontoelectron mass ratio, \(m_p/m_e\), the protontoelectron temperature ratio, \(T_p/T_e\), the ratio of some characteristic velocity, such as the electron Alfv’en speed, to the speed of light, \(v_{A_e}/c\), and the plasma beta of either the protons or the electrons, \(\beta\). It is often convenient with this normalized system to use the combination of the ratio of the electron Alfv’en speed to the speed of light and the plasma beta to derive the temperature in normalized units, like so,
assuming the plasma is quasineutral and thus, \(n_0 = 1.0\) for both the protons and electrons. The proton beta and proton temperature then follow from the specified protontoelectron temperature ratio. It is recommended that the user initialize Maxwellian distribution functions using this derived temperature, so as to avoid the ambiguity of the user’s definition of the thermal velocity,
Whether the user ultimately elects to use \(v_{th_s} = \sqrt{2 T_s/m_s}\) or \(v_{th_s} = \sqrt{T_s/m_s}\) is of no consequence to the initialization of the simulation, and likely only to manifest in the user’s specification of the velocity space extents.
From normalized to physical units in Vlasov and multifluid simulations¶
Very often we setup a problem in terms of nondimensional units. An example of one such nondimensional units is given in Normalized Units notes. Here, I describe one approach to “denormalize” the setup and create a set of dimensional values that are perhaps easier to interpret than the nondimensional units. This denormalization is not unique, expressing the fact that plasma (and fluid) equations are essentially scalefree.
As an example, consider the following fragment of a Gkeyll input file, describing a typical problem setup.
local Constants = require "Lib.Constants"
 Universal constants.
epsilon0 = 1.0  permittivity of free space.
mu0 = 1.0  permeability of free space.
lightSpeed = 1.0/math.sqrt(mu0*epsilon0)  speed of light.
 User inputs.
massRatio = 1836.153
tau = 1.0  Ratio of ion to electron temperature.
B0 = 0.25  Driven magnetic field amplitude.
beta = 0.01  Total plasma beta.
ionMass = massRatio  Ion mass in simulation.
elcMass = 1.0  Electron mass in simulation.
ionCharge = 1.0  Ion charge in simulation.
elcCharge = 1.0  Electron charge in simulation.
n = 0.01  Plasma density, same for ions and electrons.
Te = beta*(B0^2)/(2.0*mu0*(1.0+tau))  Electron temperature.
Ti = tau*Te  Ion temperature.
 Derived parameters.
vtElc = math.sqrt(2.0*Te/elcMass)
vtIon = math.sqrt(2.0*Ti/ionMass)
 cyclotron frequency
omegaCe = ionCharge*B0/elcMass
 gyro radius
rhoe = vtElc/omegaCe
Lx = 200.0*math.pi*rhoe
nuElc = 0.1*w0*(math.pi^2)  Electron collision frequency.
nuIon = nuElc/math.sqrt(ionMass*(tau^3))  Ion collision frequency.
wpe = math.sqrt(n*elcCharge^2/(elcMass*epsilon0))  plasma frequency
What do these numbers mean? For example, in the above we set \(n = 0.01\). Obviously, this does not mean that there are \(1/100\) particles per meter! If it did, then if our cellsize was say 1/10 of a meter, then we would have hardly any particles in a single cell. In fact, we could no longer treat the problem with a Vlasov equation but would need to use a Nbody description instead. So how to interpret these numbers?
First, recall that these are nondimensional quantities. Hence, these numbers do not have any units and can’t be interpreted as if they have units. To do a meaningful interpretation, we must denormalize the numbers by picking some reference values. There are several reasonable choices one can make. For example, one can choose a reference frequency. Given the speed of light in SI units we can then determine the reference length scale. Then, using SI units values for electron mass, fundamental charge and other quantities, allows us to completely denormalize all values. Of course, other reasonable choices are possible too. For example, we can select a reference length and then use the speed of light to determine the reference time (frequency) unit.
In the following fragment, we select our plasmafrequency as \(\omega_{pe} = 10^{10}\) \(s^{1}\), and the compute a reference frequency and, using speed of light, a lengthscale:
wpePhys = 1e10  Plasmafrequency in [1/s]
WDIM = wpePhys/wpe  Reference frequency [1/s]
CDIM = Constants.SPEED_OF_LIGHT  Reference speed
LDIM = CDIM/WDIM  Reference length scale
Now, all other quantities can be computed:
nPhys = wpePhys^2*Constants.ELECTRON_MASS*Constants.EPSILON0/Constants.ELEMENTARY_CHARGE^2
lambdaDPhys = vtElc/wpe*LDIM
OmegaCePhys = omegaCe*WDIM
B0Phys = OmegaCePhys*Constants.ELECTRON_MASS/Constants.ELEMENTARY_CHARGE
vAlfPhys = B0Phys/math.sqrt(Constants.MU0*Constants.ELECTRON_MASS*massRatio*nPhys)
Note we use physical SI unit values of electron mass, elementary charge, \(\epsilon_0\) and \(\mu_0\). With this, the various physical values are:
Number density 3.14208e+16 [#/m^3]
Electron thermal speed 5.29963e+06 [m/s]
Ion thermal speed 123678 [m/s]
Debye length 0.000529963 [m]
Electron gyroradius 0.000211985 [m]
Domain length 0.133194 [m]
Plasma parameter 4.67686e+06 [#]
B0 0.142141 [T]
vAlf/c 0.0583426
These numbers appear perfectly reasonable. For example, the plasma parameter, i.e. the number of particles inside a Debye sphere, is computed as \(n \lambda_D^3 = 4.7\times 10^{6}\), showing that the plasma approximation is perfectly valid.
Of course, other choices of the initial plasmafrequency (or another choice of a particular physical parameter like the domain size or numberdensity) would give a different set of values. However, of course, independent of the choice, the physics remains unchanged as long as all physical dimensions are scaled consistently. (Which is of course the virtue of the nondimensional units in the first place).
The eigensystem of the Maxwell equations with extension to perfectly hyperbolic Maxwell equations¶
Eigensystem of Maxwell equations¶
In this document I list the eigensystem of the Maxwell equations. Maxwell’s equations consist of the curl equations
along with the divergence relations
Here, \(\mathbf{E}\) is the electric field, \(\mathbf{B}\) is the magnetic flux density, \(\epsilon_0\), \(\mu_0\) are permittivity and permeability of free space, and \(\mathbf{J}\) and \(\varrho_c\) are specified currents and charges respectively. The speed of light is determined from \(c=1/(\mu_0\epsilon_0)^{1/2}\).
These are linear equations and hence the eigensytem is independent of the value of the electromagnetic fields. In 1D Maxwell equations can be written as, ignoring sources,
The eigenvalues of this system are \(\{c,c,c,c,0,0\}\). The right eigenvectors of the flux Jacobian are given by the columns of the matrix
The left eigenvectors are the rows of the matrix
Eigensystem of Perfectly Hyperbolic Maxwell equations¶
The perfectly hyperbolic Maxwell equations are a modification of the Maxwell equations that take into account the divergence relations. The modified equations explicitly “clean” divergence errors and are a hyperbolic generalization of the Hodge project method commonly used in electromagnetism to correct for charge conservation errors. See [munz_2000], [munz_2000b], [munz_2000c] for details.
These equations are written as
Here, \(\psi\) and \(\psi\) are correction potentials for the electric and magnetic field respectively and \(\chi\) and \(\gamma\) are dimensionless factors that control the speed at which the errors are propagated.
In 1D these equations can be written as, ignoring sources,
The eigenvalues of this system are \(\{c\gamma, c\gamma, c\chi, c\chi, c, c, c, c\}\). The right eigenvectors of the flux Jacobian are given by the columns of the matrix
The left eigenvectors are the rows of the matrix
[munz_2000]  C.D Munz, P. Omnes, R. Schneider and E. Sonnendruer and U. Voss, “Divergence Correction Techniques for Maxwell Solvers Based n a Hyperbolic Model”, Journal of Computational Physics, 161, 484511, 2000. 
[munz_2000b]  C.D Munz, P. Omnes, and R. Schneider, “A threedimensional finitevolume solver for the Maxwell equations with divergence cleaning on unstructured meshes”, Computer Physics Communications, 130, 83117, 2000. 
[munz_2000c]  C.D Munz and U. Voss, “A FiniteVolume Method for the Maxwell Equations in the Time Domain”, SIAM Journal of Scientific Computing, 22, 449475, 2000. 
The eigensystem of the Euler equations¶
In this document I list the eigensystem of the Euler equations valid for a general equation of state. The formulas are taken from [Kulikovskii2001], Chapter 3, section 3.1. The Euler equations can be written in conservative form as
where
is the total energy and \(\varepsilon\) is the internal energy of the fluid and \(q^2=u^2 + v^2 + w^2\). The pressure is given by an equation of state (EOS) \(p=p(\varepsilon, \rho)\). For an ideal gas the EOS is \(p = (\gamma1)\rho \varepsilon\).
The eigenvalues are \(\{uc, u, u, u, u+c\}\). The right eigenvectors of the flux Jacobian are given by the columns of the matrix
here
is the enthalpy and the sound speed respectively. Also,
Note that for ideal gas EOS we have
and \(b=\gamma1\). Hence, in this case the term \(hc^2/b\) in (2) is just \(q^2/2\). The left eigenvectors are the rows of the matrix
where
which, for an ideal gas EOS reduces to \(q^2/2\).
Now consider the problem of splitting a jump vector \(\Delta \equiv [\delta_0,\delta_1,\delta_2,\delta_3,\delta_4]^T\) into coefficients neeeded in computing the Riemann problem. The coefficients are given by \(L\Delta\). For an ideal gas law EOS, after some algebra we can show that an efficient way to compute these are
References¶
[Kulikovskii2001]  Andrei G. Kulikoviskii and Nikolai V. Pogorelov and Andrei Yu. Semenov, Mathematical Aspects of Numerical Solutions of Hyperbolic Systems, Chapman and Hall/CRC, 2001. 
The eigensystem of the tenmoment equations¶
In this document I list the eigensystem of the tenmoment equations. These equations are derived by taking moments of the Boltzmann equation and truncating the resulting infinite series of equations by assuming the heat flux tensor vanishes. In nonconservative form these equations are
In these equations square brackets around indices represent the minimal sum over permutations of free indices within the bracket needed to yield completely symmetric tensors. Note that there is one such system of equations for each species in the plasma. Here, \(q\) is the species charge, \(m\) is the species mass, \(n\) is the number density, \(u_j\) is the velocity, \(P_{ij}\) the pressure tensor and \(\mathbf{E}\) and \(\mathbf{B}\) are the electric and magnetic field respectively. Also \(\partial_t \equiv \partial /\partial t\) and \(\partial_i \equiv \partial /\partial x_i\).
To determine the eigensystem of the homogeneous part of this system we first write, in onedimension, the lefthand side of of these equations in the form
where \(\mathbf{v}\) is the vector of primitive variables and \(\mathbf{A}\) is the quasilinear coefficient matrix [1]. For the tenmoment system we have
where \(\rho \equiv mn\) and
The eigensystem of this matrix needs to be determined. It is easiest to use a computer algebra system for this. I prefer the open source package Maxima for this. The righteigenvectors returned by Maxima need to massaged a little bit to bring them into a clean form. The results are described below.
The eigenvalues of the system are given by
To maintain hyperbolicity we must hence have \(\rho>0\) and \(P_{11}>0\). In multiple dimensions, in general, the diagonal elements of the pressure tensor must be positive. When \(P_{11}=0\) the system reduces to the cold fluid equations which is known to be rank deficient and hence not hyperbolic as usually understood [2]. Also notice that the eigenvalues do not include the usual fluid soundspeed \(c_s=\sqrt{5p/3\rho}\) but instead have two different propagation speeds \(c_1=\sqrt{P_{11}/\rho}\) and \(c_2=\sqrt{3P_{11}/\rho}\). This is because the (neutral) tenmoment system does not go to the correct limit of Euler equations in the absence of collisions. In fact, it is collisions that drive the pressure tensor to isotropy. These collision terms should also be included in the plasma tenmoment system. In this case, however, the situation is complicated due to the presence of multiple species of very different masses which leads to interspecies collision terms that need to be computed carefully. For a twospecies plasma, for example, see the paper by Green [Green1973] in which the relations for relaxation of momentum and energy are used to derive a simplified collision integral for use in the Boltzmann equation.
The right eigenvectors (column vectors) are given below.
and
and
We can now compute the left eigenvectors (row vectors) by inverting the matrix with right eigenvectors stored as columns. This ensures the normalization \(\mathbf{l}^p \mathbf{r}^k = \delta^{pk}\), where the \(\mathbf{l}^p\) are the left eigenvectors. On performing the inversion we have
and
and
The eigensystem of the equations written in conservative form¶
In the wavepropagation scheme the quasilinear equations can be updated. However, the resulting solution will not be conservative. This actually might not be a problem for the tenmoment system as the system (as written) can not be put into a homogeneous conservation law form anyway. However, most often for numerical simulations the eigensystem of the conservation form of the homogeneous system is needed. This eigensystem is related to the eigensystem of the quasilinear form derived above. To see this consider a conservation law
where \(\mathbf{f} = \mathbf{f}(\mathbf{q})\) is a flux function. Now consider an invertible transformation \(\mathbf{q} = \varphi(\mathbf{v})\). This transforms the conservation law to
where \(\varphi'\) is the Jacobian matrix of the transformation and \(D\mathbf{f} \equiv \partial \mathbf{f}/\partial \mathbf{q}\) is the flux Jacobian. Comparing this to the quasilinear form we see that the quasilinear matrix is related to the flux Jacobian by
This clearly shows that the eigenvalues of the flux Jacobian are the same as those of the quasilinear matrix while the right and left eigenvectors can be computed using \(\varphi' \mathbf{r}^p\) and \(\mathbf{l}^p(\varphi')^{1}\) respectively.
For the tenmoment system the required transformation is
For this transformation we have
The inverse of the transformation Jacobian is
References¶
[Green1973]  John M. Greene. Improved BhatnagarGrossKrook model of electronion collisions. The Physics of Fluids, 16(11):20222023, 1973. 
[1]  There is no standard name for this matrix. I choose to call it the quasilinear coefficient matrix instead of the incorrect term “primitive flux Jacobian”. 
[2]  For hyperbolicity the quasilinear matrix must posses real eigenvalues and a complete set of linearly independent right eigenvectors. For the cold fluid system we only have a single eigenvalue (the fluid velocity) and a single eigenvector. This can lead to generalized solutions like delta shocks. 
Handling twofluid fivemoment and tenmoment source terms¶
The twofluid system treats a plasma as a mixture of electron and ion fluids, coupled via the electromagnetic field and collisions. In the ideal twofluid system collisions and heatflux are neglected, leading to a closed system of coupled PDEs: one set of fluid equations for each of the fluids, and Maxwell equations for the electromagnetic field. Nonneutral effects, electron interia as well as displacement currents are retained. Further, the fluid pressures can be treated as either a scalar (fivemoment model) or a symmetric \(3\times 3\) tensor (tenmoment model), or a combination.
While solving the twofluid system there are two distinct solves: the hyperbolic update, which can be performed for each equation system separately and the source update, which couples the fluids and the fields together. In this note I only focus on the source updates, for both the five as well as the tenmoment equations. When written in nonconservation law form, the only sources in the system are the Lorentz force in the momentum equation, the plasma currents in the electric field equation and, in the tenmoment model, the rotation of the pressure tensor due to the magnetic field. The latter equation is uncoupled from the source update for the momentum and the electric field, and can be treated separately.
Fivemoment source updates¶
The equations for the fivemoment source updates are
where, for the plasma species \(s\), \(n_s\) is the number density, \(\mathbf{v}_s\) is the velocity, \(q_s\) and \(m_s\) are the charge and mass respectively. Further, \(\mathbf{E}\) is the electric field, and \(\epsilon_0\) is permittivity of free space. In these equations the magnetic field and number density are constants, as there are no source terms for these quantities.
It is more convenient to work in terms of the plasma current \(\mathbf{J}_s \equiv q_s n_s \mathbf{v}_s\), which leads to the coupled system
where \(\mathbf{\Omega}_s \equiv q_s\mathbf{B}/m_s\) and \(\omega_s \equiv \sqrt{q_s^2 n_s/\epsilon_0 m_s}\) is the species plasma frequency. This is a system of linear, constantcoefficient ODEs for the \(3s+3\) unknowns \(\mathbf{J}_s\) and \(\mathbf{E}\).
Implicit solution¶
To solve the system of ODEs we replace the timederivatives with timecentered differences. This leads to the discrete equations
where \(\mathbf{J}_s^{n+1/2} = (\mathbf{J}_s^{n+1}+\mathbf{J}_s^{n})/2\) and \(\mathbf{E}^{n+1/2} = (\mathbf{E}^{n+1}+\mathbf{E}^n)/2\). This is a system of \(3p+3\) system of linear equations for the \(3p+3\) unknowns \(\mathbf{J}_s^{n+1/2}\), \(s=1,\ldots,p\) and \(\mathbf{E}_s^{n+1/2}\) and can be solved with any linear algebra routine.
Tenmoment source updates¶
The tenmoment equations have identical sources for currents and electric field. For these terms the same implicit algorithm can be used. In addition, there are source terms in the pressure equation. In nonconservative form these can be written as the linear system of equations
This system can be updated using a similar timecentered implicit method as use for the currents and the electric field. Note that unlike the sources for the currents and the electric field, the pressure source terms are uncoupled from the other fluid and electromagnetic terms.