Computational Methods for Kinetic Processes in Plasma Physics

Computational Methods for Kinetic Processes in Plasma Physics Ken Nishikawa National Space Science & Technology Center/UAH September 28, 2011 1/39

Context Electromagnetic code new_bbeps 1. f, a magnetized electromagnetic code 1. understanding the code 2. calculates 3 d phase space distribution, velocity moments, and entropy (PSDIST 13) add section to use this subroutine in newbbeps 1. f 3. contour subroutine is in nullgks 2. f

Kernel code: new_bbeps 1. f use use 1. 2. 3. 4. 5. init 1 d bpush 1 d dpush 1 d fft 1 d field 1 d diag 1 d init 1 d. mod bpush 1 d. mod dpush 1 d. mod fft 1 d. mod field 1 d. mod diag 1 d. mod init 1 mod. f bpush 1 mod. f dpush 1 mod. f fft 1 mod. f field 1 mod. f diag 1 mod. f nullgks. f init 1 lib. f (fort 77) bpush 1 lib. f dpush 1 lib. f fft 1 lib. f field 1 lib. f diag 1 lib. f understand the structure of the code understand the essential subroutines understand init 1 lib. f and init 1 mod. f understand diag 1 lib. f and diag 1 mod. f to implement new diagnostics calculates 3 d phase space distribution (PSDIST 1)

!-----------------------------------! * * * periodic 1 d electromagnetic particle simulation kernel code * * ! this is a simple 1 d skeleton particle-in-cell code designed for ! exploring new computer architectures. it contains the critical pieces ! needed for depositing charge and current, advancing particles, and ! solving the fields. the code moves only electrons, with periodic ! electromagnetic forces obtained by solving maxwell's equation with ! fast fourier transforms. ! written by viktor k. decyk, ucla ! Fortran 90 for Macintosh G 3 ! copyright 1999, regents of the university of california ! update: july 14, 2011 ! modules below are included (ken) program bbeps 1 use init 1 d initialize particle position and velocity use bpush 1 d use dpush 1 d use fft 1 d use field 1 d use diag 1 d produce plots

implicit none ! idimp = dimension of phase space = 4 ! nmv = number of segments in v for velocity distribution integer : : idimp = 4, nmv = 40, ipbc = 1 ! default unit numbers integer : : iuin = 8, iuot = 18, iudm = 19, iup = 11, iua = 15 integer : : iue = 26 integer : : np, npx 1, nxh, nxeh, nx 1 integer : : nloop, itime 0, ntime, ltime, isign, irc integer : : it, itw real : : zero = 0. 0, time = 0. 0, tloop = 0. 0 real : : tpush = 0. 0, tdpost = 0. 0, tdjpost = 0. 0, tsort = 0. 0 real : : tfft = 0. 0, totpush = 0. 0, ts = 0. 0 real : : qbme, affp, dth, qi 0, omt, q 2 m 0, wp 0 real : : we, wf, wm, wef, wke, ws real, dimension(: , : ), pointer : : part real, dimension(: ), pointer : : qe, fxe real, dimension(: , : ), pointer : : cu, amu, fxyze, byze complex, dimension(: , : ), pointer : : eyz, byz complex, dimension(: ), pointer : : ffc, ffe

integer, dimension(: ), pointer : : mixup complex, dimension(: ), pointer : : sct real, dimension(: ), pointer : : pt integer, dimension(: ), pointer : : ip, npic real, dimension(: ), pointer : : sfield complex, dimension(: ), pointer : : pott real, dimension(: , : ), pointer : : vfield complex, dimension(: , : ), pointer : : vpott, vpotr real, dimension(: , : ), pointer : : fv, fvm real, dimension(: , : ), pointer : : wt character(len=10) : : cdrun character(len=32) : : fname 991 format (' T = ', i 7) 992 format (' field, kinetic, total energies = ', 3 e 14. 7) 993 format (' electric(l, t), magnetic energies = ', 3 e 14. 7) ! read namelist iuin = get_funit(iuin) open(unit=iuin, file='input 1', form='formatted', status='old') read (iuin, input 1) read parameters in input 1 ! override input data

idcode = 2 psolve = 1 ndim = 2 ! create string from idrun write (cdrun, '(i 10)') idrun cdrun = adjustl(cdrun) ! text output file iuot = get_funit(iuot) fname = 'output 1. '//cdrun open(unit=iuot, file=trim(fname), form='formatted', status='replace') ! np = total number of electrons in simulation np = npx + npxb; npx 1 = npx + 1 npx: ambient, npxb: beam nx = 2**indx; nxh = nx/2 nxe = nx + 4 ! ax =. 866025 if (inorder==LINEAR) then ! ax =. 912871 nxe = nx + 2 endif nxeh = nxe/2

! dimension for index and sorting arrays nx 1 = nx + 1 ! nloop = number of time steps in simulation nloop = tend/dt +. 0001 ! part(1, n) = position x of particle n ! part(2, n) = velocity vx of particle n initialized by subroutine DISTR 1 H ! part(3, n) = velocity vy of particle n ! part(4, n) = velocity vz of particle n allocate(part(idimp, np)) set the dimension ! in real space, qe(j) = charge density at grid point j ! in real space, fxe(j) = longitudinal force/charge at grid point j, ! that is, fxe is the convolution of the longitudinal electric field ! over the particle shape allocate(qe(nxe), fxe(nxe)) ! in real space, fxyze(i, j) = i component of force/charge at grid (j) ! that is, fxyze are the convolutions of the electric field ! over the particle shape

allocate(fxyze(3, nxe)) ! cu(i, j) = i component of current at grid (j). ! byze(i, j) = i component of magnetic field at grid (j). ! byze is the convolution of the magnetic field over the particle shape allocate(cu(2, nxe), byze(2, nxe)) ! in fourier space, exyz = transverse electric field ! in fourier space, bxyz = magnetic field allocate(eyz(2, nxeh), byz(2, nxeh)) ! ffc = form factor array for poisson solver allocate(ffc(nxh)) ! mixup = array of bit reversed addresses for fft ! sct = sine/cosine table for fft allocate(mixup(nxh), sct(nxh)) ! ! open graphics device irc = open_graphs(nplot) enable to use GKS graphic package ! initialize timer call wtimer(time, ltime, -1)

! initialize constants itime 0 = 0 itime = itime 0 ntime = itime + itime 0 qbme = qme affp = real(nx)/real(np) omt = sqrt(omy*omy + omz*omz) q 2 m 0 = qbme*qme*real(np)/real(nx) dth = 0. 0 ! debug ! dth =. 5*dt ! end debug wp 0 = q 2 m 0*affp plasma frequency ! set initial time t 0 = dt*real(itime 0) ! set default diagnostic file names if (ntp > 0) fpname = 'potk 1. '//cdrun if (nta > 0) faname = 'vpotk 1. '//cdrun if (nte > 0) fename = 'vpotrk 1. '//cdrun

! energy diagnostics if (ntw > 0) then allocate(wt((nloop-1)/ntw-(itime 0/ntw)+1, 7)) itw = 0 endif ! initialize electromagnetic fields byze = 0. 0 if (omt > 0. 0) then call baddext(byze, omy, omz, nx, inorder) call cguard(byze, nx, inorder) endif byz = cmplx(0. 0, 0. 0) eyz = cmplx(0. 0, 0. 0) cu = 0. 0 ! prepare fft tables call fft_init(mixup, sct, indx) ! calculate form factors call pois_init(ffc, ax, affp, nx)

! initialize density profile and velocity distribution ! background electrons if (npx > 0) call distr(part, 1, npx, vty, vtz, vx 0, vy 0, vz 0, npx, nx, &ipbc) ! beam electrons if (npxb > 0) call distr(part, npx 1, npxb, vtdx, vtdy, vtdz, vdx, vdy, vdz& &, npxb, nx, ipbc) Remark! in init 1 mod. f distr is assined by `distr 1 h’ ! distr => idistrh 1 initializes x and vx, vy, vz co-ordinates for ! magnetized 1 -2/2 d codes. ! calls DISTR 1 H (see the next page) ! initialize charge density to background qi 0 = -qme/affp ! fix guiding centers for electrons if (omt > 0. 0) call distr(part, byze, np, qbme, nx, ipbc, inorder) ! retard electron positions to deposit current ! call retard(part, np, dth, nx) ! sorting arrays if (sortime > 0) allocate(pt(np), ip(np), npic(nx 1)) ! initialize diagnostics ! diagnostic metafile iudm = get_funit(iudm)

! velocity diagnostic if (ntv > 0) then allocate(fv(2*nmv+2, 3), fvm(3, 3)) fv(1, : ) = 8. *max(vtx, vty, vtz) endif ! potential diagnostic if (ntp > 0) then allocate(sfield(nxe)) if (modesxp > nxh) modesxp = nxh allocate(pott(modesxp)) ! open output file if (nprec==0) then nprec = -1; iup = get_funit(iup) call bfopen(pott, modesxp, iup, nprec, trim(fpname)) endif ! vector potential or electromagnetic diagnostic if ((nta > 0). or. (nte > 0)) then allocate(vfield(2, nxe)) vfield = 0. 0 endif

! vector potential diagnostic if (nta > 0) then if (modesxa > nxh) modesxa = nxh allocate(vpott(2, modesxa)) ! open output file if (narec==0) then narec = -1; iua = get_funit(iua) call bfopen(vpott, modesxa, iua, narec, trim(faname)) endif ! electromagnetic diagnostic if (nte > 0) then if (modesxe > nxh) modesxe = nxh allocate(vpotr(2, modesxe)) ! open output file if (nerec==0) then nerec = -1; iue = get_funit(iue) call bfopen(vpotr, modesxe, iue, nerec, trim(fename)) endif

! record time call wtimer(time, ltime) write (iuot, *) 'initialization wall clock time = ', time, 'sec' ! ! * * * start main iteration loop * * * ! 500 if (nloop <= ntime) go to 2000 write (iuot, 991) ntime ! ! prepare electromagnetic diagnostic if (nte > 0) then it = ntime/nte if (ntime==nte*it) vfield = cu endif ! initialize current density to background call sguard(cu, zero, nx, inorder) ! deposit electron current call djpost(part, cu, np, qme, dth, tdjpost, nx, ipbc, inorder, djopt)

! initialize charge density to background call sguard(qe, qi 0, nx, inorder) ! deposit electron charge call dpost(part, qe, np, qme, tdpost, inorder, dopt) ! add guard cells for current call aguard(cu, nx, inorder) ! add guard cells call aguard(qe, nx, inorder) ! velocity diagnostic if (ntv > 0) then for example ntv=5 in input 1. whistler it = ntime/ntv if (ntime==ntv*it) then ! calculate paticle distribution function and moments call vdist(part, fvm, np, nmv) ! display velocity distributions call displayfv(fv, fvm, ' ELECTRON', ntime, nmv, 2, irc) if (irc==1) go to 2000 endif

! phase velocity diagnostic if (ntv > 0) then for example ntv=5 in input 1. whistler it = ntime/ntv if (ntime==ntv*it) then ! calculate 3 d phase space distribution abd moments call PSDIST 13(part, fpsm, nxb, idimp, nmv, nmvf) How can we plot contour plots for phase space?

! transform charge to fourier space isign = -1 call fft(qe, isign, mixup, sct, tfft, indx, inorder) ! potential diagnostic if (ntp > 0) then it = ntime/ntp if (ntime==ntp*it) then ! calculate potential in fourier space isign = 1 call pois(qe, sfield, isign, ffc, ws, nx, inorder) ! store selected fourier modes call gtmodes(sfield, pott, nx, modesxp, inorder) ! write diagnostic output call writebf(pott, modesxp, iup, nprec, order=LINEAR) ! transform potential to real space call fft(sfield, isign, mixup, sct, tfft, indx, inorder) call cguard(sfield, nx, inorder)

! display potential call displays(sfield, ' POTENTIAL', ntime, 999, 0, nx, irc, inorder& &) if (irc==1) go to 2000 endif ! transform current to fourier space isign = -1 call fft(cu, isign, mixup, sct, tfft, indx, inorder) ! electromagnetic diagnostic if (nte > 0) then it = ntime/nte if (ntime==nte*it) then ! calculate averaged radiative vector potential vfield = 0. 5*(vfield + cu) call avrpot(vfield, byz, ffc, ci, nx, inorder) ! store selected fourier modes call gtmodes(vfield, vpotr, nx, modesxe, inorder) ! write diagnostic output call writebf(vpotr, modesxe, iue, nerec, order=LINEAR)

! transform radiative vector potential to real space isign = 1 call fft(vfield, isign, mixup, sct, tfft, indx, inorder) call cguard(vfield, nx, inorder) ! display radiative vector potential call displayv(vfield, ' RADIATIVE VPOTENTIAL', ntime, 999, 0, 1, & &nx, irc, inorder) if (irc==1) go to 2000 endif ! calculate electromagnetic fields in fourier space if (ntime==0) then ! calculate initial darwin magnetic field call ibpois(cu, byz, ffc, ci, wm, nx, inorder) wf = 0. ! calculate initial darwin electric field allocate(amu(2, nxe), ffe(nxeh))

! deposit momentum flux call sguard(amu, zero, nx, inorder) call dmjpost(part, amu, np, qme, ts, inorder, djopt) call aguard(amu, nx, inorder) ! solve for darwin electric field isign = -1 call fft(amu, isign, mixup, sct, tfft, indx, inorder) ! byze = cu call dcuperp(byze, amu, nx, inorder) call epois_init(ffe, ax, affp, wp 0, ci, nx) call iepois(byze, eyz, ffe, ci, wf, nx, inorder) deallocate(amu, ffe) dth =. 5*dt ! calculate electromagnetic fields else call maxwel(eyz, byz, cu, ffc, ci, dt, wf, wm, nx, inorder) endif

! vector potential diagnostic if (nta > 0) then it = ntime/nta if (ntime==nta*it) then ! calculate vector potential in fourier space call avpot(byz, vfield, nx, inorder) ! store selected fourier modes call gtmodes(vfield, vpott, nx, modesxa, inorder) ! write diagnostic output call writebf(vpott, modesxa, iua, narec, order=LINEAR) ! transform vector potential to real space isign = 1 call fft(vfield, isign, mixup, sct, tfft, indx, inorder) call cguard(vfield, nx, inorder) ! display vector potential call displayv(vfield, ' VECTOR POTENTIAL', ntime, 999, 0, 1, nx, & &irc, inorder) if (irc==1) go to 2000 endif

! calculate longitudinal electric field in fourier space isign = -1 call pois(qe, fxe, isign, ffc, we, nx, inorder) ! add longitudinal and transverse electric fields call emfield(fxyze, fxe, eyz, ffc, nx, inorder) ! copy magnetic field call emfield(byze, byz, ffc, nx, inorder) ! transform force/charge to real space isign = 1 call fft(fxyze, isign, mixup, sct, tfft, indx, inorder) call cguard(fxyze, nx, inorder) ! transform magnetic field to real space isign = 1 call fft(byze, isign, mixup, sct, tfft, indx, inorder) ! add external magnetic field if (omt > 0. 0) call baddext(byze, omy, omz, nx, inorder) call cguard(byze, nx, inorder) ! push particles wke = 0. call push 3(part, fxyze, byze, omx, np, qbme, dth, wke, tpush, nx, ipbc, & &inorder, popt)

! sort electrons if (sortime > 0) then if (mod(ntime, sortime)==0) then call sortp(part, pt, ip, npic, tsort, inorder) endif ! energy diagnostic if (ntw > 0) then it = itime/ntw if (itime==ntw*it) then wef = we + wf + wm ws = wef + wke write (iuot, 992) wef, wke, ws write (iuot, 993) we, wf, wm itw = itw + 1 wt(itw, : ) = (/wef, wke, zero, ws, we, wf, wm/) endif itime = itime + 1 ntime = itime + itime 0

call wtimer(tloop, ltime) time = time + tloop go to 500 2000 continue ! ! * * * end main iteration loop * * * ! ! energy diagnostic if (ntw > 0) then ts = t 0 + dt*real(ntw)*(itime 0 -(itime 0/ntw)*ntw) call reset_graphs call displayw(wt, ts, dt*real(ntw), itw, irc) endif ! accumulate timings write (iuot, *) 'electromagnetic code bbeps 1' write (iuot, *) 'main wall clock time = ', time, 'sec' totpush = tpush + tdpost + tdjpost write (iuot, *) 'electron push time = ', tpush, 'sec' write (iuot, *) 'electron charge deposit time = ', tdpost, 'sec'

write (iuot, *) 'electron current deposit time = ', tdjpost, 'sec' write (iuot, *) 'total electron push time = ', totpush, 'sec' write (iuot, *) 'electron sort time = ', tsort totpush = totpush + tsort write (iuot, *) 'total electron time = ', totpush, 'sec' write (iuot, *) 'total fft time=', tfft, 'sec' time = time - (totpush + tfft) write (iuot, *) 'other time=', time, 'sec' ! write final diagnostic metafile fname = 'diag 1. '//cdrun open(unit=iudm, file=trim(fname), form='formatted', status='replace') ! potential diagnostics if (ntp > 0) then nprec = nprec - 1 ceng = affp write (iudm, pot 1 d) endif

! vector potential diagnostics if (nta > 0) then narec = narec - 1 ceng = affp write (iudm, vpot 1 d) endif ! electromagnetic diagnostics if (nte > 0) then nerec = nerec - 1 ceng = affp write (iudm, em 1 d) endif ! write out input file write (iudm, input 1) write (iuot, *) ' * * * q. e. d. * * *' close(unit=iudm) close(unit=iuot) ! close graphics device call close_graphs stop end program bbeps 1

diag 1 lib. f c-----------------------------------c 1 d PIC library for diagnostics c diag 1 lib. f contains diagnostic procedures: c VDIST 1 calculates 1 component velocity distribution, velocity moments, c and entropy. c VDIST 13 calculates 3 component velocity distribution, velocity c moments, and entropy. c PSDIST 1 calculates 1 d phase space distribution, velocity moments, and c entropy. c PSDIST 13 calculates 3 d phase space distribution, velocity moments, and c entropy. attention c FWRITE 1 writes real binary data to file c FCWRITE 1 writes complex binary data to file. c FREAD 1 writes real binary data to file. c FCREAD 1 writes complex binary data to file. c written by viktor k. decyk, ucla c copyright 1994, regents of the university of california c update: july 7, 2010 c------------------------------------

subroutine VDIST 13(part, fvm, idimp, nmv, nmvf) c for 1 -2/2 d code, this subroutine calculates 3 d velocity distribution, c velocity moments, and entropy c input: all except fvm, output: fv, fvm c part(2, n) = velocity vx of particle n c part(3, n) = velocity vy of particle n c part(4, n) = velocity vz of particle n c fv = distribution function, number of particles in each velocity range c maximum velocity (used for scaling) is contained in first element fv. c vdrift for i-th dimension is contained in fvm(1, i) c vth for i-th dimension is contained in fvm(2, i) c entropy for i-th dimension is contained in fvm(3, i), defined to be: c s/k = -sum(f(v)/np)*log(f(v)/(np*delta_v)). Assumes that distribution c is uniform in space and distributions in each dimension are c independent. c idimp = size of phase space = 2 c np = number of particles

c nmv = number of segments in v for velocity distribution c the number of velocity bins used is 2*nmv + 1, nmvf >= 2*nmv+2 (nmv=40) implicit none integer idimp, nmv, nmvf real part, fvm dimension part(idimp, np), fv(nmvf, 3), fvm(3, 3) c local data double precision sumvx, sumvy, sumvz, sumvx 2, sumvy 2, sumvz 2, anp real anmv, svx, svy, svz integer j, nvx, nvy, nvz anmv = real(nmv) c maximum velocity (used for scaling) is contained in first element fv(1, 1), - svx = anmv/fv(1, 1) svy = anmv/fv(1, 2) svz = anmv/fv(1, 3) c zero out distribution do 10 j = 2, nmvf fv(j, 1) = 0. 0 fv(j, 2) = 0. 0 fv(j, 3) = 0. 0 10 continue

c count particles in each velocity region anmv = anmv + 2. 5 sumvx = 0. 0 d 0 sumvy = 0. 0 d 0 sumvz = 0. 0 d 0 sumvx 2 = 0. 0 d 0 sumvy 2 = 0. 0 d 0 sumvz 2 = 0. 0 d 0 do 20 j = 1, np nvx = part(2, j)*svx + anmv sumvx = sumvx + part(2, j) sumvx 2 = sumvx 2 + part(2, j)**2 nvy = part(3, j)*svy + anmv sumvy = sumvy + part(3, j) sumvy 2 = sumvy 2 + part(3, j)**2 nvz = part(4, j)*svz + anmv sumvz = sumvz + part(4, j) sumvz 2 = sumvz 2 + part(4, j)**2 if ((nvx. ge. 2). and. (nvx. le. nmvf)) fv(nvx, 1) = fv(nvx, 1) + 1. 0 if ((nvy. ge. 2). and. (nvy. le. nmvf)) fv(nvy, 2) = fv(nvy, 2) + 1. 0 if ((nvz. ge. 2). and. (nvz. le. nmvf)) fv(nvz, 3) = fv(nvz, 3) + 1. 0 20 continue

c calculate velocity moments anp = 1. 0 d 0/dble(np) sumvx = sumvx*anp fvm(1, 1) = sumvx fvm(2, 1) = dsqrt(sumvx 2*anp - sumvx**2) sumvy = sumvy*anp fvm(1, 2) = sumvy fvm(2, 2) = dsqrt(sumvy 2*anp - sumvy**2) sumvz = sumvz*anp fvm(1, 3) = sumvz fvm(2, 3) = dsqrt(sumvz 2*anp - sumvz**2) c calculate entropy sumvx = 0. 0 d 0 sumvy = 0. 0 d 0 sumvz = 0. 0 d 0 sumvx 2 = 0. 0 d 0 sumvy 2 = 0. 0 d 0 sumvz 2 = 0. 0 d 0

do 30 j = 2, nmvf if (fv(j, 1). gt. 0. 0) then sumvx = sumvx + fv(j, 1) sumvx 2 = sumvx 2 + fv(j, 1)*dlog(dble(fv(j, 1)*svx)) endif if (fv(j, 2). gt. 0. 0) then sumvy = sumvy + fv(j, 2) sumvy 2 = sumvy 2 + fv(j, 2)*dlog(dble(fv(j, 2)*svy)) endif if (fv(j, 3). gt. 0. 0) then sumvz = sumvz + fv(j, 3) sumvz 2 = sumvz 2 + fv(j, 3)*dlog(dble(fv(j, 3)*svz)) endif 30 continue if (sumvx. gt. 0. 0 d 0) sumvx = -sumvx 2/sumvx + dlog(sumvx) if (sumvy. gt. 0. 0 d 0) sumvy = -sumvy 2/sumvy + dlog(sumvy) if (sumvz. gt. 0. 0 d 0) sumvz = -sumvz 2/sumvz + dlog(sumvz) fvm(3, 1) = sumvx fvm(3, 2) = sumvy fvm(3, 3) = sumvz return end

subroutine PSDIST 13(part, fpsm, nxb, idimp, nmv, nmvf) c for 1 -2/2 d code, this subroutine calculates 3 d phase space attention c distribution, velocity moments, and entropy c input: all except fpsm and psm, output: fps, fpsm, psm c part(1, n) = position x of particle n c part(2, n) = velocity vx of particle n c part(3, n) = velocity vy of particle n c part(4, n) = velocity vz of particle n c fps = distribution function, number of particles in each velocity c and spatial range c maximum velocity (used for scaling) is contained in first element fps. c vdrift for i-th dimension is contained in fpsm(1, i) c vth for i-th dimension is contained in fpsm(2, i) c entropy for i-th dimension is contained in fpsm(3, i), defined to be: c s/k = -sum(f(v)/np)*log(f(v)/(np*delta_v)). Assumes that c distributions in each dimension are independent. c psm = double precision scratch array c nx = system length in x direction c nxb = number of spatial regions for distribution c idimp = size of phase space = 2

c np = number of particles c the number of velocity bins used is 2*nmv + 1, nmvf >= 2*nmv+2 implicit none integer nx, nxb, idimp, nmv, nmvf real part, fpsm dimension part(idimp, np), fps(nmvf, 3, nxb), fpsm(3, 3, nxb) double precision psm dimension psm(2, 3, nxb) c local data double precision sumvx, sumvy, sumvz, sumvx 2, sumvy 2, sumvz 2, anp real anmv, svx, svy, svz, at 1 integer j, k, n, nvx, nvy, nvz anmv = real(nmv) c zero out distribution do 20 k = 1, nxb do 10 j = 2, nmvf fps(j, 1, k) = 0. 0 fps(j, 2, k) = 0. 0 fps(j, 3, k) = 0. 0 10 continue

fpsm(1, 1, k) = anmv/fps(1, 1, k) fpsm(1, 2, k) = anmv/fps(1, 2, k) fpsm(1, 3, k) = anmv/fps(1, 3, k) psm(1, 1, k) = 0. 0 d 0 psm(2, 1, k) = 0. 0 d 0 psm(1, 2, k) = 0. 0 d 0 psm(2, 2, k) = 0. 0 d 0 psm(1, 3, k) = 0. 0 d 0 psm(2, 3, k) = 0. 0 d 0 20 continue c count particles in each velocity and spatial region at 1 = real(nxb)/real(nx) anmv = anmv + 2. 5 do 30 j = 1, np n = part(1, j)*at 1 n=n+1 if ((n. ge. 1). and. (n. le. nxb)) then nvx = part(2, j)*fpsm(1, 1, n) + anmv psm(1, 1, n) = psm(1, 1, n) + part(2, j) psm(2, 1, n) = psm(2, 1, n) + part(2, j)**2 nvy = part(3, j)*fpsm(1, 2, n) + anmv

psm(1, 2, n) = psm(1, 2, n) + part(3, j) psm(2, 2, n) = psm(2, 2, n) + part(3, j)**2 nvz = part(4, j)*fpsm(1, 3, n) + anmv psm(1, 3, n) = psm(1, 3, n) + part(4, j) psm(2, 3, n) = psm(2, 3, n) + part(4, j)**2 if ((nvx. ge. 2). and. (nvx. le. nmvf)) fps(nvx, 1, n) = fps(nvx, 1, n) + 1 1. 0 if ((nvy. ge. 2). and. (nvy. le. nmvf)) fps(nvy, 2, n) = fps(nvy, 2, n) + 1 1. 0 if ((nvz. ge. 2). and. (nvz. le. nmvf)) fps(nvz, 3, n) = fps(nvz, 3, n) + 1 1. 0 endif 30 continue c calculate velocity moments anp = 1. 0 d 0/dble(np) do 50 k = 1, nxb svx = fpsm(1, 1, k) svy = fpsm(1, 2, k) svz = fpsm(1, 3, k) sumvx = psm(1, 1, k)*anp

fpsm(1, 1, k) = sumvx fpsm(2, 1, k) = dsqrt(psm(2, 1, k)*anp - sumvx**2) sumvy = psm(1, 2, k)*anp fpsm(1, 2, k) = sumvy fpsm(2, 2, k) = dsqrt(psm(2, 2, k)*anp - sumvy**2) sumvz = psm(1, 3, k)*anp fpsm(1, 3, k) = sumvz fpsm(2, 3, k) = dsqrt(psm(2, 3, k)*anp - sumvz**2) c calculate entropy sumvx = 0. 0 d 0 sumvy = 0. 0 d 0 sumvz = 0. 0 d 0 sumvx 2 = 0. 0 d 0 sumvy 2 = 0. 0 d 0 sumvz 2 = 0. 0 d 0 do 40 j = 2, nmvf if (fps(j, 1, k). gt. 0. 0) then sumvx = sumvx + fps(j, 1, k) sumvx 2 = sumvx 2 + fps(j, 1, k)*dlog(dble(fps(j, 1, k)*svx)) endif

if (fps(j, 2, k). gt. 0. 0) then sumvy = sumvy + fps(j, 2, k) sumvy 2 = sumvy 2 + fps(j, 2, k)*dlog(dble(fps(j, 2, k)*svy)) endif if (fps(j, 3, k). gt. 0. 0) then sumvz = sumvz + fps(j, 3, k) sumvz 2 = sumvz 2 + fps(j, 3, k)*dlog(dble(fps(j, 3, k)*svz)) endif 40 continue if (sumvx. gt. 0. 0 d 0) sumvx = -sumvx 2/sumvx + dlog(sumvx) if (sumvy. gt. 0. 0 d 0) sumvy = -sumvy 2/sumvy + dlog(sumvy) if (sumvz. gt. 0. 0 d 0) sumvz = -sumvz 2/sumvz + dlog(sumvz) fpsm(3, 1, k) = sumvx fpsm(3, 2, k) = sumvy fpsm(3, 3, k) = sumvz 50 continue return end

!-----------------------------------! module diag 1 d ! ! Fortran 90 interface to 1 d PIC Fortran 77 library diag 1 lib. f ! diag 1 mod. f contains diagnostic procedures: ! defines module diag 1 d ! wtimer performs wall clock timing. ! get_funit returns an unconnected fortran unit number. ! bfopen => bfopen 1 opens binary file for real 1 d scalar data. ! bfopen => bfvopen 1 opens binary file for real 1 d vector data ! bfopen => bfcopen 1 opens binary file for complex 1 d scalar data. ! bfopen => bfvcopen 1 opens binary file for complex 1 d vector data ! writebf => ifwrite 1 writes real 1 d scalar data to a binary file. ! calls FWRITE 1 ! writebf => ifvwrite 1 writes real 1 d vector data to a binary file. ! call FWRITE 1 ! writebf => ifcwrite 1 writes complex 1 d scalar data to a binary file. ! calls FCWRITE 1

! writebf => ifvcwrite 1 writes complex 1 d vector data to a binary file. ! call FCWRITE 1 ! readbf => ifread 1 reads real 1 d scalar data from a binary file. ! calls FREAD 1 ! readbf => ifvread 1 reads real 1 d vector data from a binary file. ! calls FREAD 1 ! readbf => ifcread 1 reads complex 1 d scalar data from a binary file. ! calls FCREAD 1 ! readbf => ifvcread 1 reads complex 1 d vector data from a binary file. ! calls FCREAD 1 ! vdist => ivdist 1 calculates 1 or 3 component velocity distribution, ! velocity moments, and entropy. ! calls VDIST 13 or VDIST 1 ! psdist => ipsdist 1 calculates 1 or 3 component phase space ! distribution, velocity moments, and entropy ! calls PSDIST 13 or PSDIST 1 attention ! written by viktor k. decyk, ucla ! copyright 2000, regents of the university of california ! update: july 18, 2011

! use globals, only: LINEAR, QUADRATIC use init 1 d, only: idrun, indx, ntp, nta, nte, psolve, tend, dt, & &ndim, omx, omy, omz, ci, ntd, ntj, t 0, ceng, pot 1 d, modesxp, & &nprec, fpname, vpot 1 d, modesxa, narec, faname, em 1 d, modesxe, &nerec, fename, den 1 d, modesxd, ndrec, fdname, vcur 1 d, modesxj, &njrec, fjname ! implicit none private public : : LINEAR, QUADRATIC public : : idrun, indx, ntp, nta, nte, psolve public : : tend, dt, ndim, omx, omy, omz, ci, ntd, ntj public : : t 0, ceng public : : pot 1 d, modesxp, nprec, fpname public : : vpot 1 d, modesxa, narec, faname public : : em 1 d, modesxe, nerec, fename public : : den 1 d, modesxd, ndrec, fdname public : : vcur 1 d, modesxj, njrec, fjname & &

! public : : wtimer, get_funit public : : bfopen, writebf, readbf public : : writef public : : vdist, psdist public : : open_graphs, close_graphs, reset_graphs public : : displays, displayv, displayfv, displayw ! ! npl = (0, 1) = display is (off, on) integer, save : : npl = 1 ! ! define interface to original Fortran 77 procedures ! interface subroutine FCWRITE 1(f, nxv, iunit, nrec, lrec, name) implicit none integer : : nx, nxv, iunit, nrec, lrec ! complex, dimension(*) : : f complex : : f character(len=*) : : name end subroutine end interface

interface subroutine VDIST 13(part, fvm, idimp, nmv, nmvf) integer : : idimp, nmv, nmvf real, dimension(idimp, np) : : part real, dimension(nmvf, 3) : : fv real, dimension(3, 3) : : fvm end subroutine end interface subroutine PSDIST 13(part, fpsm, nxb, idimp, nmv, nmvf) integer : : nx, nxb, idimp, nmv, nmvf attention real, dimension(idimp, np) : : part real, dimension(nmvf, 3, nxb) : : fps real, dimension(3, 3, nxb) : : fpsm double precision, dimension(2, 3, nxb) : : psm end subroutine end interface

interface subroutine DISPR(f, label, xmin, xmax, isc, ist, mks, nxv, ngs, chr, c& &hrs, irc) implicit none integer : : isc, ist, mks, nxv, ngs, irc real : : xmin, xmax character(len=*) : : label, chr character(len=*), dimension(ngs) : : chrs ! real, dimension(*) : : f real : : f end subroutine end interface ! ! define generic interfaces to Fortran 90 library ! interface bfopen module procedure bfopen 1 module procedure bfvopen 1 module procedure bfcopen 1 module procedure bfvcopen 1 end interface

! interface writebf module procedure ifwrite 1 module procedure ifvwrite 1 module procedure ifcwrite 1 module procedure ifvcwrite 1 end interface ! interface readbf module procedure ifread 1 module procedure ifvread 1 module procedure ifcread 1 module procedure ifvcread 1 end interface ! interface writef module procedure ifwrite 0 end interface

! interface vdist module procedure ivdist 1 end interface ! interface psdist module procedure ipsdist 1 end interface ! interface displays module procedure idscaler 1 end interface ! interface displayv module procedure idvector 1 end interface ! interface displayfv module procedure idisplayfv 1 end interface

! interface displayw module procedure idisplayw 1 end interface ! subroutine ivdist 1(part, fvm, np, nmv) attention ! calculates 1 d velocity distribution, velocity moments, and entropy integer : : np, nmv real, dimension(: , : ), pointer : : part real, dimension(: , : ), pointer : : fv, fvm ! local data integer : : idimp, nmvf, idimv idimp = size(part, 1) nmvf = size(fv, 1); idimv = size(fv, 2) if (idimv==1) then call VDIST 1(part, fvm, idimp, nmv, nmvf) else if (idimv==3) then call VDIST 13(part, fvm, idimp, nmv, nmvf) endif end subroutine ivdist 1

! subroutine ipsdist 1(part, fpsm, nx, np, nmv) attention ! calculates 1 d phase space distribution, velocity moments, and entropy integer : : nx, np, nmv real, dimension(: , : ), pointer : : part real, dimension(: , : ), pointer : : fps, fpsm ! local data double precision, dimension(2, size(fps, 2), size(fps, 3)) : : psm integer : : idimp, nmvf, idimv, nxb idimp = size(part, 1) nmvf = size(fps, 1); idimv = size(fps, 2); nxb = size(fps, 3) if (idimv==1) then call PSDIST 1(part, fpsm, nxb, idimp, nmv, nmvf) else if (idimv==3) then call PSDIST 13(part, fpsm, nxb, idimp, nmv, nmvf) endif end subroutine ipsdist 1

! subroutine idisplayfv 1(fv, fvm, label, itime, nmv, idt, irc) ! displays velocity distribution functions ! fv = velocity distribution ! fvm = velocity moments ! itime = current time step ! nmv = number of velocity intervals ! idt = (1, 2, 3) = display (individual, composite, both) functions idt=2 ! irc = return code (0 = normal return) implicit none integer : : itime, nmv, idt, irc real, dimension(: , : ), pointer : : fv, fvm character(len=*) : : label ! isc = 999 = use default display scale ! ist = 1 = display positive values only ! mks = 0 = cycle through line styles integer : : isc = 999, ist = 1, mks = 0 integer : : i, nmvf, nmv 2, idimv real : : vmax, vmin

character(len=12) : : c character(len=2) : : cs character(len=54) : : lbl character(len=45) : : chr character(len=10), dimension(3) : : chrs 91 format(', T =', i 7) 92 format(' VD =', f 9. 6, ' VTH =', f 9. 5) 93 format(' VTX =', f 9. 5) 94 format(' VTX =', f 9. 5, ' VTY =', f 9. 5, ' VTZ =', f 9. 5) ! chrs = short array of characters to label individual line samples data chrs /' VX ', ' VY ', ' VZ '/ if (npl==0) return idimv = size(fv, 2) if ((idimv /= 1). and. (idimv /= 3)) return nmvf = size(fv, 1) nmv 2 = 2*nmv + 1 write (c, 91) itime

! each velocity distributions on its own plot if (idt /= 2) then do i = 1, idimv cs = trim(adjustl(chrs(i))) lbl = trim(label)//' VELOCITY DISTR VS '//cs//c write (chr, 92) fvm(1, i), fvm(2, i) vmax = fv(1, i) vmin = -vmax call DISPR(fv(2, i), lbl, vmin, vmax, isc, ist, mks, nmv 2, nmvf, 1, chr& &, chrs(i), irc) if (irc > 127) then npl = irc - 128 if (npl==0) call CLRSCRN irc = 0 return endif if (irc==1) return enddo endif

! all velocity distributions on common plot if (idt /= 1) then lbl = trim(label)//' VELOCITY DISTRS VS '//'V'//c if (idimv==1) then write (chr, 93) fvm(2, 1) vmax = fv(1, 1) vmin = -vmax else if (idimv==3) then write (chr, 94) fvm(2, 1), fvm(2, 2), fvm(2, 3) vmax = max(fv(1, 1), fv(1, 2), fv(1, 3)) vmin = -vmax endif call DISPR(fv(2, 1), lbl, vmin, vmax, isc, ist, mks, nmv 2, nmvf, idimv& &, chrs, irc) if (irc > 127) then npl = irc - 128 if (npl==0) call CLRSCRN irc = 0 endif end subroutine idisplayfv 1

c Null general 1 d gks graphics library c written by viktor k. decyk, ucla c copyright 1994, regents of the university of california c update: march 21, 2009 c-----------------------------------subroutine GROPEN c this subroutine opens gks and activates standard workstation c colors and maximum size of display surface are also set c if a string (keyboard) device is available, it is initialized c if a locator (mouse) device is available, it is initialized return end c-----------------------------------subroutine GROPEN 0(iecho) c this subroutine opens gks and activates workstation c colors and maximum size of display surface are also set c if a string (keyboard) device is available, it is initialized c if a locator (mouse) device is available, it is initialized return end

c-----------------------------------subroutine GRCLOSE c this subroutine deactivates workstation and closes gks return end c-----------------------------------subroutine SETNPLT(nplt, irc) c this subroutine resets the maximum number of plots per page c if requested number is negative, it is set to default (=1) c the current plot location is also reset to the next available location c if next available location is iplot = 0 and the old location was not c 0, then the workstation is updated and a return code can be generated. irc = 0 return end

in nullgks 1. f c-----------------------------------subroutine DISPR(f, label, xmin, xmax, isc, ist, mks, nxv, ngs, chrs 1, irc) c this subroutine plots ngs subarrays of the array f, on a common graph, c each plot with nx points, versus a linear function in x, c where xmin < xmax. c depending on the number of colors in the display device, each subarray c is plotted with a different color, given in order by: c blue, red, yellow, cyan, magenta, green, and foreground. c after all the colors are cycled through, then different line styles c are cycled through if mks=0, in order: solid, dash, dot, dash-dot, c or different marker types if mks>0: dot, plus, star, circle, cross. c multiple plots per page can be displayed by dividing the screen into c n x n subregions, where n*n is the next largest integer >= nplot c the location (ix, iy) of a plot in the subregions is determined by c the parameter iplot = ix + iy*n c f = array containing subarrays to be plotted c label = long character string label for plot c xmin/xmax = range of x values in plot

c isc = power of 2 scale of y coordinate for plot c ist = flag for choosing positive and/or negative values c the plots have a common scale in y given by ymax and ymin. c if ist = 0, then ymax = 2**isc and ymin = -2**isc. c if ist = 1, then ymax = 2**isc and ymin = 0. c if ist = -1, then ymax = 0 and ymin = -2**isc. c if ist = 2, then ymin = fmin, ymax = fmin + 2**ir, c where fmin/fmax are the function minimum/maximum, c and ir = power of 2 scale for (fmax - fmin) c if abs(isc) < 116, then the isc value passed is used for scale. c if abs(isc) > 116, then the program finds the minimum value of isc c which will contain the plots, determined by the absolute value of f c mks = flag to determine whether lines or markers are used, c mks=0 means cycle through lines styles, mks > 0 means cycle through c marker styles, using the value of mks as the initial marker, c mks < 0 means draw the first subarray with a line, then subsequent c subarrays with markers, using the value of abs(mks) as the initial c marker. c nx = number of points plotted in each subarray c nxv = first dimension of array f, nxv >= nx

c ngs = second dimension of array f, number of subarrays to be plotted c chr = additional long character string comment for plot c chrs = array of ngs short character labels used by subroutine tickd c to label individual line or marker samples c irc = return code (0 = normal return) c nxbs = length of scratch variable for plotting character*(*) label, chr character*(*) chrs(ngs) dimension f(nxv, ngs) irc = 0 return end Find the plot subroutine for color contours (see vspectrum. f)?

Summary in new_bbeps 1. f ! calculate paticle distribution function and moments call vdist(part, fvm, np, nmv) ! display velocity distributions call displayfv(fv, fvm, ' ELECTRON', ntime, nmv, 2, irc) in diag 1 lib. f subroutine VDIST 13(part, fvm, idimp, nmv, nmvf) in nullgks. f subroutine DISPR(f, label, xmin, xmax, isc, ist, mks, nxv, ngs, chrs 1, irc) in diag 1 mod. f ! subroutine idisplayfv 1(fv, fvm, label, itime, nmv, idt, irc)

! displays velocity distribution functions ! fv = velocity distribution ! fvm = velocity moments ! itime = current time step ! nmv = number of velocity intervals ! idt = (1, 2, 3) = display (individual, composite, both) functions (idt=2) ! irc = return code (0 = normal return) implicit none integer : : itime, nmv, idt, irc real, dimension(: , : ), pointer : : fv, fvm character(len=*) : : label ! isc = 999 = use default display scale ! ist = 1 = display positive values only ! mks = 0 = cycle through line styles integer : : isc = 999, ist = 1, mks = 0 integer : : i, nmvf, nmv 2, idimv real : : vmax, vmin character(len=12) : : c character(len=2) : : cs

character(len=54) : : lbl character(len=45) : : chr character(len=10), dimension(3) : : chrs 91 format(', T =', i 7) 92 format(' VD =', f 9. 6, ' VTH =', f 9. 5) 93 format(' VTX =', f 9. 5) 94 format(' VTX =', f 9. 5, ' VTY =', f 9. 5, ' VTZ =', f 9. 5) ! chrs = short array of characters to label individual line samples data chrs /' VX ', ' VY ', ' VZ '/ if (npl==0) return idimv = size(fv, 2) if ((idimv /= 1). and. (idimv /= 3)) return nmvf = size(fv, 1) nmv 2 = 2*nmv + 1 write (c, 91) itime

! each velocity distributions on its own plot if (idt /= 2) then do i = 1, idimv cs = trim(adjustl(chrs(i))) lbl = trim(label)//' VELOCITY DISTR VS '//cs//c write (chr, 92) fvm(1, i), fvm(2, i) vmax = fv(1, i) vmin = -vmax call DISPR(fv(2, i), lbl, vmin, vmax, isc, ist, mks, nmv 2, nmvf, 1, chr& &, chrs(i), irc) if (irc > 127) then npl = irc - 128 if (npl==0) call CLRSCRN irc = 0 return endif if (irc==1) return enddo endif

! all velocity distributions on common plot if (idt /= 1) then lbl = trim(label)//' VELOCITY DISTRS VS '//'V'//c if (idimv==1) then write (chr, 93) fvm(2, 1) vmax = fv(1, 1) vmin = -vmax else if (idimv==3) then write (chr, 94) fvm(2, 1), fvm(2, 2), fvm(2, 3) vmax = max(fv(1, 1), fv(1, 2), fv(1, 3)) vmin = -vmax endif call DISPR(fv(2, 1), lbl, vmin, vmax, isc, ist, mks, nmv 2, nmvf, idimv& &, chrs, irc) if (irc > 127) then npl = irc - 128 if (npl==0) call CLRSCRN irc = 0 endif end subroutine idisplayfv 1

! interface subroutine DISPR(f, label, xmin, xmax, isc, ist, mks, nxv, ngs, chr, c& &hrs, irc)