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 