Skeleton 2-1/2D Darwin Particle-in-Cell (PIC) codes
by Viktor K. Decyk
copyright 2007-2015, regents of the university of california

This program contains sample codes for illustrating the basic structure
of a 2-1/2D Darwin 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 PIC codes.  They are also intended as benchmark
reference codes to aid in developing new codes and in evaluating new
computer architectures.

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 DModels.pdf.

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 near-field electromagnetic interaction, obtained
by solving the Darwin subset of Maxwell's equation, where the transverse
displacement current in neglected in Ampere's law.  A spectral method
using Fast Fourier Transforms (FFTs) is used to solve the scalar and
vector 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.
Additional details can be found in the companion presentations,
Darwin.pdf and Darwin2.pdf.

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, a current
derivative, and charge deposit, add guard cell procedures, a vector and
scalar FFT, transverse current procedures, Poisson solvers, vector FFTs,
copy guard cell procedures, a particle push, and a particle sorting
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:
   GJPOST2L (cgjpost2l): deposit current density
   GPOST2L (cgpost2l): deposit charge density
   GDJPOST2L (cgdjpost2l): deposit acceleration density and momentum
                           flux
   GDCJPOST2L (cgdcjpost2l): deposit electron current and acceleration
                             density and momentum flux
   ASCFGUARD2L (cascfguard2l): add acceleration density and scaled
                               transverse electric field
   AGUARD2L (caguard2l): add charge density guard cells
   ACGUARD2L (cacguard2l): add current density guard cells
   ACGUARD2L (cacguard2l): add acceleration density guard cells
   AMCGUARD2L (camcguard2l): add momemtum flux guard cells

Field solve section:
   WFFT2RX (cwfft2rx): FFT charge density to fourier space
   POIS23 (cpois23): calculate smoothed longitudinal electric field in
                     fourier space.
   WFFT2R3 (cwfft2r3): FFT smoothed longitudinal electric field to real
                       space
   WFFT2R3 (cwfft2r3): FFT current density to fourier space
   CUPERP2 (ccuperp2): take transverse part of current
   BBPOIS23 (cbpois23): calculate magnetic field in fourier space
   WFFT2R3 (cwfft2r3): FFT smoothed magnetic field to real space
   BADDEXT2 (cbaddext2): add constant to magnetic field
   WFFT2R3 (cwfft2r3): FFT acceleration density to fourier space
   WFFT2RN (cwfft2rn): FFT momentum flux to fourier space
   ADCUPERP23 (cadcuperp23): take transverse part of time derivative of
                             current from momentum flux and acceleration
                             density
   EPOIS23 (cepois23): calculate convective part of transverse electric
                       field
   WFFT2R3 (cwfft2r3): FFT smoothed transverse electric field to real
                       space

Particle Push section:
   BGUARD2L (cbguard2l): fill in guard cells for smoothed longitudinal
                         electric field
   BGUARD2L (cbguard2l): fill in guard cells for smoothed magnetic field
   BGUARD2L (cbguard2l): fill in guard cells for smoothed transverse
                         electric field
   ADDVRFIELD2 (caddvrfield2): add longitudinal and transverse electric
                               fields
   GBPUSH23L (cgbpush23l): update particle co-ordinates with smoothed
                           electric and magnetic fields:
                           x(t)->x(t+dt/2); v(t-dt/2)->v(t+dt/2)
   DSORTP2YL (cdsortp2yl) : sort particles by cell

The inputs to the code are the grid parameters indx, indy, the particle
number parameters npx, npy, the time parameters tend, dt, and the
velocity parameters vtx, vty, vtz, vx0, vy0, vz0, the inverse speed of
light ci, magnetic field electron cyclotron frequencies omx, omy, omz,
the sorting parameter sortime, and number of corrections in darwin
iteration ndc

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.
   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.
   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
sortime = number of time steps between electron sorting.
   This is used to improve cache performance.  sortime=0 to suppress.
omx/omy/omz = magnetic field electron cyclotron frequency in x/y/z
ndc = number of corrections in darwin iteration
   typical values are 1 or 2.

The major program files contained here include:
dpic2.f90    Fortran90 main program 
dpic2.c      C main program
dpush2.f     Fortran77 procedure library
dpush2_h.f90 Fortran90 procedure interface (header) library
dpush2.c     C procedure library
dpush2.h     C 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 .c
and .h adhere to the C99 standard.

The makefile is setup to use gcc and gfortran with Linux.  A version for
Mac OS X is also in the Makefile, but is commented out.  

Two executables can be created, fdpic2 for Fortran and cdpic2 for C.

To compile program, execute:

Make program_name

where program_name is either: fdpic2 or cdpic2, or execute:

make

to create both programs.

To execute, type the name of the executable:

./program_name

where program_name is either fdpic2 or cdpic2.

The file output contains the results produced for the default parameters.

The Fortran version can be compiled to run with double precision by
changing the Makefile (typically by setting the compiler options flags
-r8).

The library dpush2.c contains wrapper functions to allow the C library
to be called from Fortran. The library dpush2_f.c contains wrapper
functions to allow the Fortran library to be called from C.