Skeleton 3D Electromagnetic OpenMP/Vector Particle-in-Cell (PIC) codes
by Viktor K. Decyk
copyright 2016, regents of the university of california

This program contains sample codes for illustrating the basic structure
of a 3D Electromagnetic OpenMP/Vector Particle-in-Cell (PIC) code, in
both Fortran and C. The codes have no diagnosics except for initial and
final energies.  Their primary purpose is to provide example codes for
physical science students learning about OpenMP/Vector PIC codes.  They
are also intended as benchmark reference codes to aid in developing new
codes and in evaluating new computer architectures.  An OpenMP/serial
version of this code with the same structure (mbpic3) also exists, and
can be compared to this code in order to understand the parallel
algorithms.

PIC codes are widely used in plasma physics.  They model plasmas as
particles which interact self-consistently via the electromagnetic
fields they themselves produce.  PIC codes generally have three
important procedures in the main iteration loop.  The first is the
deposit, where some particle quantity, such as a charge, is accumulated
on a grid via interpolation to produce a source density.  The second
important procedure is the field solver, which solves Maxwell’s equation
or a subset to obtain the electric and/or magnetic fields from the
source densities.  Finally, once the fields are obtained, the particle
forces are found by interpolation from the grid, and the particle
co-ordinates are updated, using Newton’s second law and the Lorentz
force.  The particle processing parts dominate over the field solving
parts in a typical PIC application.

More details about PIC codes can be found in the texts by C. K. Birdsall
and A. B. Langdon, Plasma Physics via Computer Simulation, 1985,
R. W. Hockney and J. W. Eastwood, Computer Simulation Using Particles,
1981, and John M. Dawson, "Particle simulation of plasmas", Rev. Mod.
Phys. 55, 403 (1983).  Details about the mathematical equations and
units used in this code is given in the companion article,
"Description of Electromagnetic Spectral Code from the UPIC Framework"
by Viktor K. Decyk, UCLA, in the file EMModels.pdf.

Details abut OpenMP can be found in the book by Rohit Chandra, Leonardo
Dagum, Dave Kohr, Dror Maydan, Jeff McDonald, and Ramesh Menon, Parallel
Programming in OpenMP, Morgan Kaufmann, 2001.

Vectorization today is largely based on exploiting Single Instruction
Multiple Data (SIMD) processors that can perform the same instruction
on multiple data items.  The one we are focussing on here is the Intel
PHI Coprocessor Knight's Corner (KNC) Many Integrated Core (MIC).  Two
approaches are used for vectorization.  The simplest is to use
vectorization tools that will analyze and compile Fortran and C/C++
codes to run on the SIMD processors.  Alternatively, Intel also provides
vector intrinsincs, which are assembly-coded functions that allow direct
access to vector assembly language instructions in C/C++.  Compiler
vectorization is normally implemented from loops.  Vector instrinics are
implemented as a vector language. Further details about Intel compiler
vectorization can be found at:
https://software.intel.com/en-us/intel-vectorization-tools
Further details about Intel vector intrinsics can be found at:
https://software.intel.com/sites/products/documentation/doclib/iss/2013/
compiler/cpp-lin/index.htm#GUID-28F49E4A-615D-4373-A36E-C8A92E913325.htm
and at https://software.intel.com/sites/landingpage/IntrinsicsGuide/.

No warranty for proper operation of this software is given or implied.
Software or information may be copied, distributed, and used at own
risk; it may not be distributed without this notice included verbatim
with each file.  If use of these codes results in a publication, an
acknowledgement is requested.

The code here uses the complete electromagnetic interaction, obtained by
solving Maxwell's equation.  A spectral method using Fast Fourier
Transforms (FFTs) is used to solve the Maxwell and Poisson equations.  A
real to complex FFT is used, and the data in Fourier space is stored in
a packed format, where the input and output sizes are the same.  The
boundary conditions are periodic, only electron species are included,
and linear interpolation is used.

For parallelization, the code uses a tiling (or blocking) technique.
Space is divided into small 3D tiles (with typically 8x8x8 grid points
in a tile), and particles are organized so that their co-ordinates in x,
y, and z lie within the same tile and are stored together in memory.
Assigning different threads to different tiles avoids data collisions
(where 2 threads try to write to the same memory location at the same
time).  The size of the tiles should be smaller than the L1 cache and
the number of tiles should be equal to or greater than the number of
processing cores.  There are 4 major procedures which make use of the
tiled data structure, the charge and current deposits, the particle
push, and the particle reordering.  In the deposit procedures, each
thread first deposits to a small local density array the size of a tile
plus guard cells.  After all the particles in a tile have been processed
the small density array is added to the global density arrays.  The
particle push is similar, where each thread first copies the global
field arrays to small local field arrays which are then used for field
interpolation.  The particle reordering procedure (VPPORDERF3LT) moves
particles which are leaving one tile into the appropriate location in
another tile.  This is done by first copying outgoing particles into an
ordered particle buffer.  Then particles are copied from the buffer into
the new locations in the particle array.  Performing this in two steps
allows one to avoid data collisions when running in parallel.  Further
information about this tiling parallel algorithm used can be found in
the companion presentation OpenMP-PIC.pdf and in the article:
V. K. Decyk and T. V. Singh, "Particle-in-Cell Algorithms for Emerging
Computer Architectures," Computer Physics Communications, 185, 708,
2014, available at http://dx.doi.org/10.1016/j.cpc.2013.10.013.

Both vectorization approaches are illustrated in this Skeleton code.
One version uses the KNC MIC vector instrinsics.  Two libraries,
avx512lib3.c and kncmbpush3.c, implement the functionality of the
OpenMP/serial mbpush3.c file in the mbpic3 directory.  This is very low
level programming and difficult to understand, but gives the best
performance.  The second version made use of compiler directives.  Two
libraries, vmbpush3.f and vmbpush3.c, reimplemented the OpenMP/serial
mbpush3.f and mbpush3.c files in such a way that the compiler could
automatically vectorize the loops.  Details of about this process for
this code are described in the file mVectorPIC3.pdf.  A parameter kvec
in the main codes selects which version will run.

The default particle push does not calculate the list of particles
leaving the tiles.  This was done because the code was faster.  There
is, however, a version of the push (VGBPPUSHF3LT) which does calculate
the list.  This version of the push requires a different reordering
procedure (VPPORDERF3LT) which does not calculates the list.

Particles are initialized with a uniform distribution in space and a
gaussian distribution in velocity space.  This describes a plasma in
thermal equilibrium.  The inner loop contains a current and charge
deposit, add guard cell procedures, a vector and scalar FFT, a
transverse current procedure, Maxwell and Poisson solvers, vector FFTs,
copy guard cell procedures, a particle push, and a particle reordering
procedure.  The final energy and timings are printed.  A sample output
file for the default input parameters is included in the file output.

In more detail, the inner loop of the code contains the following
procedures in Fortran (C):

Deposit section:
   VGRJPPOST3LT (cvgrjppost3lt): relativistic current deposit, and
            or cknc2grjppost3lt  update position. x(t+dt/2)->x(t+dt)
   VGJPPOST3LT (cvgjppost3lt): deposit current density, and update
           or cknc2gjppost3lt  position. x(t+dt/2)->x(t+dt)
   VPPORDER3LT (cvpporder3lt) : move particles to appropriate tile
            or ckncpporder3lt
   VGPPOST3LT (cvgppost3lt): deposit charge density
          or cknc2gppost3lt
   ACGUARD3L (cacguard3l): add current density guard cells
         or ckncacguard3l
   AGUARD3L (caguard3l): add charge density guard cells
        or ckncaguard3l

Field solve section:
   WFFT3RVMX (cwfft3rvmx): FFT charge density to fourier space
          or ckncwfft3rmx
   WFFT3RVM3 (cwfft3rvm3): FFT current density to fourier space
          or ckncwfft3rm3
   MCUPERP3 (cmcuperp3): take transverse part of current
        or ckncmcuperp3
   VMIBPOIS33 (cvmibpois33): calculate initial magnetic field in fourier
           or ckncmibpois33  space
   VMMAXWEL3 (cvmmaxwel3): update transverse electromagnetic fields in
          or ckncmmaxwel3  fourier space
   VMPOIS33 (cvmpois33): calculate smoothed longitudinal electric field
         or ckncmpois33  in fourier space.
   VMEMFIELD3 (cvmemfield3): add smoothed longitudinal and transverse
           or ckncmemfield3  electric fields
   VMEMFIELD3 (cvmemfield3): copy smoothed magnetic field
           or ckncmemfield3
   WFFT3RVM3 (cwfft3rvm3): FFT smoothed electric field to real space
          or ckncwfft3rm3
   WFFT3RVM3 (cwfft3rvm3): FFT smoothed magnetic field to real space
          or ckncwfft3rm3

Particle Push section:
   CGUARD3L (ccguard3l): fill in guard cells for smoothed electric field
        or cknccguard3l
   CGUARD3L (ccguard3l): fill in guard cells for smoothed magnetic field
        or cknccguard3l
   V2GRBPPUSH3LT (cv2grbppush3lt): update relativistic particle
            or ckncgrbppush3lt    co-ordinates with smoothed electric
                                   and magnetic fields.
                                   x(t)->x(t+dt/2); v(t-dt/2)->v(t+dt/2)
   V2GBPPUSH3LT (cvgbppush3lt): update particle co-ordinates with
             or ckncgbppush3lt  smoothed electric and magnetic fields.
                                x(t)->x(t+dt/2); v(t-dt/2)->v(t+dt/2)
   VPPORDER3LT (cvpporder3lt) : move particles to appropriate tile
            or ckncpporder3lt

The inputs to the code are the grid parameters indx, indy, indz, the
particle number parameters npx, npy, npz, the time parameters tend, dt,
and the velocity paramters vtx, vty, vtz, vx0, vy0, vz0, the inverse
speed of light ci, the flag relativity.  In addition, a tile size mx,
my, mz, and overflow size xtras are defined as well as the version
selector kvec.

In more detail:
indx = exponent which determines length in x direction, nx=2**indx.
indy = exponent which determines length in y direction, ny=2**indy.
indz = exponent which determines length in z direction, nz=2**indz.
   These ensure the system lengths are a power of 2.
npx = number of electrons distributed in x direction.
npy = number of electrons distributed in y direction.
npz = number of electrons distributed in z direction.
   The total number of particles in the simulation is npx*npy.
tend = time at end of simulation, in units of plasma frequency.
dt = time interval between successive calculations.
   The number of time steps the code runs is given by tend/dt.
   dt should be less than .45*ci for the electromagnetic code.
vtx/vty/vtz = thermal velocity of electrons in x/y/z direction.
   a typical value is 1.0.
vx0/vy0/vz0 = drift velocity of electrons in x/y/z direction.
ci = reciprocal of velocity of light
relativity = (no,yes) = (0,1) = relativity is used
mx/my/mz = number of grids points in x, y, and z in each tile
   should be less than or equal to 16.
xtras = fraction of extra particles needed for particle management
kvec = (1,2) = run (autovector,KNC) version

The major program files contained here include:
vmbpic3.f90       Fortran90 main program 
vmbpic3.f03       Fortran2003 main program 
vmbpic3.c         C main program
omplib.f          Fortran77 OpenMP utility library
omplib_h.f90      Fortran90 OpenMP utility interface (header) library
omplib.c          C OpenMP utility library
omplib.h          C OpenMP utility header library
vmbpush3.f        Fortran77 procedure library
vmbpush3_h.f90    Fortran90 procedure interface (header) library
vmbpush3.c        C procedure library
vmbpush3.h        C procedure header library
avx512lib3.c      C Vector intrinsics utility library
avx512lib3.h      C Vector intrinsics utility header library
avx512flib3.c     C Vector intrinsics utility library for Fortran90
                  arrays
avx512lib3_h.f90  Fortran90 Vector intrinsics utility header library
avx512flib3_h.f90 Fortran90 Vector intrinsics utility header library
                  for Fortran90 arrays
avx512lib3_c.f03  Fortran2003 Vector intrinsics utility header library
kncmbpush3.c      C Vector intrinsics procedure library
kncmbpush3.h      C Vector intrinsics procedure header library
kncmbpush3_h.f90  Fortran90 Vector intrinsics procedure header library
kncmbpush3_c.f03  Fortran2003 Vector intrinsics procedure header library
dtimer.c          C timer function, used by both C and Fortran

Files with the suffix.f90 adhere to the Fortran 90 standard, files with
the suffix .f adhere to the Fortran77 standard, files with the suffix
.f03 adhere to the Fortran2003 standard, and files with the suffix .c
and .h adhere to the C99 standard.

The makefile is setup to use icc and ifort with Linux.  Definitions for
other compilers are in the Makefile, but are commented out. 

Three executables can be created, fvmbpic3 for Fortran90, f03vmbpic3 for
Fortran2003, and cvmbpic3 for C.  The only differences between the
Fortran90 and Fortran2003 codes is how the interface to C is handled.
Interoperability with C is part of Fortran2003, but is still unfamiliar
to many programmers.

To compile the program for execution in native mode, execute:

Make program_name

where program_name is either: fvmbpic3, f03vmbpic3 or cvmbpic3

To create all of them, execute:

make

To execute, first copy the executable to the Coprocessor:

scp program_name your_userid_id@mic0:program_name

then ssh to the Coprocessors and type the name of the executable:

ssh mic0
./program_name

where program_name is either fvmbpic3, f03vmbpic3 or cvmbpic3

By default, OpenMP will use the maximum number of processors it can
find.  If the user wants to control the number of threads, a parameter
nvp can be set in the main program.  In addition, the environment
variable OMP_NUM_THREADS may need to be set to the maximum number of 
threads expected.  Finally, if OpenMP 3.1 is available, setting the
environment variable OMP_PROC_BIND=true generally gives better
performance by preventing threads from moving between CPUs.

The file output contains the results produced for the default parameters
The file output.knc contains result for the KNC version.

The Fortran version can be compiled to run with double precision by
changing the Makefile (typically by setting the compiler options flags
-r8), but only when the runtime parameter kvec = 1.

The libraries omplib.c and vmbpush3.c contain wrapper functions to allow
the C libraries to be called from Fortran.  The libraries omplib_f.c and
vmbpush3_f.c contain wrapper functions to allow the Fortran libraries to
be called from C.
