Skeleton 3D Darwin OpenMP Particle-in-Cell (PIC) codes field diagnostics
by Viktor K. Decyk
copyright 2007-2015, regents of the university of california

This directory contains procedures that are useful for performing wave
diagnostics, such as measuring the wavenumber and frequency of plasma
waves, for the 3D Darwin OpenMP Particle-in-Cell (PIC) codes.  The
procedures are written in Fortran in the file mdfield3.f.  The file
mdfield3_f.c contains wrapper functions to allow the Fortran procedures
to be called from C.

The electric and magnetic fields used to advance particles are smoothed
by suppressing large wavenumbers in order to avoid numerical errors due
to aliasing.  The field diagnostics, however, are normally unsmoothed.
The main diagnostics are the electrostatic potential and the vector
potential in fourier space, in the Coulomb gauge.  These are calculated
by the procedures MPOTP3 and MAPOTP33, respectively.  Other diagnostics,
such as the unsmoothed longitudinal electric field and magnetic field
can be calculated from the potentials with the procedures MGRADF3 and
MCURLF3, respectively.  The procedure MEPOIS33 can be called with
isign=1 to obtain the unsmoothed transverse electric field.  The fourier
modes are normally stored to disk for later post-processing, such as by
performing a fourier transform in time to obtain the frequencies for
different wavenumbers.  It is usually not necessary to store all the
modes to disk, since the high wavenumber modes are usually unimportant.
Furthermore, the fourier modes are stored in a packed format which is
difficult to work with.  Thus procedures RDMODES3 and RDVMODES3 are also
provided to extract a subset of the fourier modes for the potential and
vector fields respectively, and copy them in an unpacked format for
easier post-processing. (The procedures WRMODES3 and WRVMODES3 are
provided to repack the scalar or vector fourier modes, if needed.)

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.

To use these procedures, one first needs to add some declarations in the
main programs mdpic3.90 and mdpic3.c in the parent directory mdpic3,
similar to the following.

For Fortran:
   integer :: modesx, modesy, modesz, modesxd, modesyd, modeszd
   complex, dimension(:,:,:), pointer :: potc, pott
   complex, dimension(:,:,:,:), pointer :: vpotc, vpott

For C:
   int modesx, modesy, modesz, modesxd, modesyd, modeszd;
   float complex *potc = NULL, *pott = NULL;
   float complex *vpotc = NULL, *vpott = NULL;
   
Then before the main iteration loop, one needs to allocate the arrays.

For Fortran:
   modesx = nx/4; modesy = ny/4; modesz = nz/4
   modesxd = modesx
   modesyd = min(2*modesy-1,ny); modeszd = min(2*modesz-1,nz)
   allocate(potc(nxeh,nye,nze),vpotc(ndim,nxeh,nye,nze))
   allocate(pott(modesxd,modesyd,modeszd))
   allocate(vpott(ndim,modesxd,modesyd,modeszd))

For C:
   modesx = nx/4; modesy = ny/4; modesz = nz/4;
   modesxd = modesx; 
   modesyd = 2*modesy - 1; modeszd = 2*modesz - 1;
   modesyd = modesyd < ny ? modesyd : ny;
   modeszd = modeszd < nz ? modeszd : nz;
   potc = (float complex *) malloc(nxeh*nye*nze
           *sizeof(float complex));
   vpotc = (float complex *) malloc(ndim*nxeh*nye*nze
            *sizeof(float complex));
   pott = (float complex *) malloc(modesxd*modesyd*modeszd
           *sizeof(float complex));
   vpott = (float complex *) malloc(ndim*modesxd*modesyd*modeszd
            *sizeof(float complex));

Finally, in the main iteration loops just before the push procedures,
one can add some diagnostics, such as:

For Fortran:
! perform potential diagnostic and unpack into array pott
! updates potc, pott, wt
   call MPOTP3(qe,potc,ffc,wt,nx,ny,nz,nxeh,nye,nze,nxh,nyh,nzh)
   call RDMODES3(potc,pott,nx,ny,nz,modesx,modesy,modesz,nxeh,nye,nze,&
  &modesxd,modesyd,modeszd)
! perform vector potential diagnostic and unpack into array vpott
! updates vpotc, vpott, wt
   call MAPOTP33(cue,vpotc,ffc,ci,wt,nx,ny,nz,nxeh,nye,nze,nxh,nyh,nzh)
   call RDVMODES3(vpotc,vpott,nx,ny,nz,modesx,modesy,modesz,ndim,nxeh,&
  &nye,nze,modesxd,modesyd,modeszd)

For C:
/* perform potential diagnostic and unpack into array pott */
/* updates potc, pott, wt                                  */
   cmpotp3((float complex *)qe,potc,ffc,&wt,nx,ny,nz,nxeh,nye,nze,nxh,
           nyh,nzh);
   crdmodes3(potc,pott,nx,ny,nz,modesx,modesy,modesz,nxeh,nye,nze,
             modesxd,modesyd,modeszd);
/* perform vector potential diagnostic and unpack into array vpott */
/* updates vpotc, vpott, wt                                        */
   cmapotp33((float complex *)cue,vpotc,ffc,ci,&wt,nx,ny,nz,nxeh,nye,
             nze,nxh,nyh,nzh);
   crdvmodes3(vpotc,vpott,nx,ny,nz,modesx,modesy,modesz,ndim,nxeh,nye,
              nze,modesxd,modesyd,modeszd);

One would have to modify the Makefile as well to include the files
mdfield3.f and mdfield3_f.c, as needed.
