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 