where filename is the name of the metafile to be plotted. Similarly, the metafiles can be viewed independently of the main simulation by typing:

After the simulation is finished, the user is then invited to run a post-processor, or quit. Running the post-processors is similar to run- ning the simulation program.

When the user is done with the simulation, he or she should sign off SUNSPOT with the following command:

It is also possible to execute the post-processors independently of the main simulation. For example, to execute the potential post- processor, type:

where filename is the name of the file containing the potential data. If you haven't given it a name, it is called pot1.d. If the output of the potential post-processor is sent to a metafile, it will be called pgraph.m, unless you give it a different name, and the user can then plot or view it if he or she wished.

To execute the debye cloud profile post-processor, type:

where sfilename is the name of the file containing the current potential data and rfilename is the name of the file containing the reference pot- tial. One of these files must have been given a name, and the other one one might be called pot1.d. If the output of the debye cloud post- processor is sent to a metafile, it will be called cgraph.m, unless you give it a different name.

To execute the spatial profile post-processor, type:

where filename is the name of the file containing the profile data. If you haven't given it a name, it is called pro1.d. If the output of the profile post-processor is sent to a metafile, it will be called fgraph.m, unless you give it a different name.

1. INTRODUCTION

The basic loop of the program is quite simple. From the known ini- tial positions of the particles, x(i), one calculates the charge den- sity q(x), at the grid locations x, by interpolation using quadratic spline functions. Schematically,

where S is the particle shape function, and the sum is over particles i. The electric field E(x) is then calculated at the grid points x by solv- ing Poisson's equation,

using an FFT algorithm. The force on each particle, at its location x(i), is found by interpolating back from the known electric field. Finally, the particle trajectories are updated by integrating Newton's second law,

using a time-centered leap-frog scheme. The cycle then repeats.

From the particle co-ordinates, a number of useful macroscopic (i.e., average) quantities, such as potential or velocity distribution func- tions, are calculated as diagnostics along the way to tell us about the behavior of the computer plasma.

In this code, all lengths are normalized to the grid space d, which will be related to a physical length later. The physical size of the simulation is then given by Lx = d*NX, where NX is the number of grid points used (NX = 2**INDX). All times are normalized to the inverse electron plasma frequency Wpe**-1, where Wpe**2 = 4*pi*n0*e**2/Me and n0 is the total density of electrons, background electrons plus beam. The velocities are normalized to the combination d*Wpe. Thus in one unit of computer time, a particle with one unit of velocity moves one grid point. In dimensionless units, the electron debye length Ld/d = Vt/d*Wpe, is equal to the dimensionless electron thermal velocity VTX, which is one of the input parameters to the simulation. The grid space d is thus related to the electron debye length through the relation d = Ld/VTX.

The physical electric field E is normalized by the expression, e*E/Me*(Wpe**2)*d. In these units, a unit electric field will increase an electron's velocity by one velocity unit in one unit of computer time. The physical potential V is similarly normalized by the expres- sion e*V/Me*(Wpe**2)*d**2. In these units one can show that the dimen- sionless potential is the same as (e*V/k*Te)*VTX**2.

Physical energies W are normalized to W/Me*(Wpe**2)*d**2. Thus one unit of computer energy is twice the kinetic energy of an electron moving with unit velocity.

The program can also be run relativistically if the parameter CI > 0. In that case, all quantities which previously referred to velocity, such as VTX, now refer to the momentum divided by the rest mass. The equations of motion are then:

2. MAIN PART OF SIMULATION

1. The charge deposition is performed as follows:

where x is the gridpoint nearest to x(i), and dx(i) = x(i) - x.

2. The electric field is solved in two steps. First a periodic solu- tion is found, which satisfies the inhomogeneous equation, but the wrong boundary conditions. To this is added a solution of Laplace's equation to give the correct boundary conditions. The periodic solution is found by first taking the Fast Fourier Trans- form (FFT) of the charge density:

where k = 2*pi*n/NX, n integer, S(k) is the Fourier Transform of the particle shape function, and the sum is over particles i. The periodic electric field Ep(x) can then be found from the ex- pression:

where the sum is over positive and negative values of n (zero ex- cluded). In one dimension the general solution of Laplace's equa- tion is just a linear function, so the final solution (allowing for non-neutral plasmas) is:

where the constant A must be determined by the boundary condi- tions.

3. The new particle co-ordinates are found from the expressions:

where the electric field at the particle's location is given by:

and x is the gridpoint nearest to x(i), and dx(i) = x(i) - x. If the program is relativistic, then the new co-ordinates are found from:

3. DIAGNOSTIC PART OF SIMULATION

1. The trajectory diagnostic stores the co-ordinates of a selected group of test charges and displays a subset of these at the end of the run. Either ions or electrons can be chosen as test charges.

2. The velocity-space diagnostic is actually three separate diagnos- tics. First it calculates and displays the velocity distribution function, in separate regions of space, if desired. This is cal- culated like the charge distribution, by accumulating number den- sity on a velocity-space grid. Second this diagnostic calculates and displays spatial profiles of density, kinetic energy, current, and/or kinetic energy flux. Third a scatter plot of particle phase space (v versus x for all particles) is displayed. This diagnostic is performed for both ions and electrons.

3. The ion density diagnostic stores and displays the ion density, which is useful for looking at ion waves.

4. The potential diagnostic stores and displays the potential, which is useful for looking at plasma waves.

5. The energy diagnostic displays the electric field energy, the ki- netic energy of each species, the net energy of the system (total field and particle energy minus any work done by time-varying boundary conditions), the momentum of each species and the total momentum of the system.

4. INPUT VARIABLES

• IBCS = 1 means the potential is specified at both boundaries.
  
• IBCS = 2 means vacuum is specified at both boundaries.
  
• IBCS = 3 means vacuum is specified at the left boundary, and the
  
• potential is specified at the right.
  
• IBCS = 4 means the potential is specified at the left boundary,
  and the electric field is specified at the right.
  
• IBCS = 5 means vacuum is specified at the left boundary, and the
  electric field is specified at the right.
  
• IBCS = 6 means the electric field is specified at both boundaries.
  
• IBCS = 7 means periodic boundary conditions.

Note that a special case occurs when IBCS = 7, and BVL is not zero. Namely, an external pump field is added to the self- consistent periodic electric field, of the form:

where K0 = 2*PI*AMODEX/NX, and E0 is determined by:

BVR: The boundary value for the field at x = NX. BVR can mean either the potential V(x=NX), the normal electric field EX(x=NX), or 4*pi times an external surface charge density SIG(x=NX), depend- ing on the boundary conditions specified by the parameter IBCS. Specifically, for

Note that a special case occurs when IBCS = 7, and BVR is not zero. Namely, an external pump field is added to the self- consistent periodic electric field, of the form:

where K0 = 2*PI*AMODEX/NX, and E0 is determined by:

TEND: The time at the end of the simulation, in units of Wpe**-1.

DT: The integration time step, which should satisfy the condition Wmax*DT<<1, where Wmax is the maximum frequency one expects. For plasma oscillations, DT = .2 is a typical value.

RMASS: If ions are moving, RMASS is the ion/electron mass ratio.

RTEMP: If ions are moving, RTEMP is the electron/ion initial temper- ature ratio. This ratio is used to fix the initial ion ther- mal velocity.

QMI: The charge on an ion, in units of e. A typical value is 1.

VDX: The drift (average) velocity for the secondary subgroup of elec- trons.

VTDX: The initial thermal velocity for the secondary subgroup of electrons.

QMB: The charge on a beam electron, in units of e. A typical value is -1.

NTI: The simulation can push ions with a larger time step than elec- trons since the ions do not respond to high frequency fields. NTI is the number of electron time steps between ion time steps. Values larger than about 10 should be avoided. (NTI = 0 means the two time step scheme is not used.)

AX: Each particle in the simulation has a finite shape, determined by the convolution of a gaussian filter,

with the interpolation function, where AX is the width of the filter in normalized units. If one uses quadratic spline inter- polation and has AX > .5, then this convolution is a wider gaus- sian with an effective width given by sqrt(.25 + AX*AX). To give reasonable suppression of grid errors while maintaining accurate physics, a typical choice is AX = .5*sqrt(3), which gives an ef- fective particle size of one grid space.

EDGE: For non-periodic boundary conditions, particles are specularly reflected when their spatial co-ordinate x(i) is less than EDGE or greater than NX - EDGE. This parameter is also used to initially distribute the particles in space.

NTR: The number of electron time steps between the writing of restart data on an external file, for use in either continuing a run at a later time, or restarting a run which has abnormally terminat- ed. (NTR = 0 means a restart file is not created.)

The meanings of the variables in the second (diagnostic) subgrouping are as follows:

NPROBT: The number of test charges whose trajectories will be stored.

NDIST: The number of test charges (out of NPROBT) whose x co- ordinates as a function of time will be displayed at the end of the simulation. The first NDIST particles out of NPROBT will be displayed.

MODES: The number of Fourier modes calculated by the potential diag- nostic.

NDP: The number of times the potential diagnostic is performed before the potential in real space is displayed. In other words, the potential is displayed every NTP*NDP electron time steps. If NDP = 0, the potential is never displayed, but is still calcula- ted and stored to disk.

where N is an integer < 5, S(x) is the particle shape function, and the sum is over particles. Thus density = P0(x), current = P1(x), energy = .5*Mass*P2(x), and energy flux = .5*Mass*P3(x). Depending on the values of JPRO and NDV, these profiles are also displayed. Finally, a phase space plot of the particle co- ordinates (v versus x) is displayed, depending on the values of JPS and NDV.

NMV: The number of segments into which the positive velocity range is divided when calculating the velocity distribution function. The positive velocity range is given by (0,VM), where VM is the power of 2 which is greater than or equal to the velocity magni- tude of the fastest particle. The negative velocity range is given by (0,-VM). Thus the size of the velocity bin is given by VM/NMV. NMV = 40 is a typical value.

NXB: The number of regions into which space is divided for the pur- pose of calculating velocity distribution functions. Velocity distributions are calculated separately in each spatial region, given by NX*IX/NXB < x < NX*(IX+1)/NXB, where IX is an integer < NXB.

NPRS: The starting index for spatial profiles. NPRO profiles are calculated, starting with the profile corresponding to index N = NPRS, and ending with the profile corresponding to index N = NPRO - NPRS + 1. The profile index N = 1 corresponds to density, N = 2 to energy, N = 3 to current, and N = 4 to energy flux.

NPRO: The number of spatial profiles to be calculated.

JPRO: Flag which determines whether the spatial profiles are ever displayed (JPRO = 1), or not (JPRO = 0). They are still calcu- lated and stored to disk, however, if NTV > 0.

JPS: Flag which determines whether phase space plots of the particles are ever displayed (JPS = 1), or not (JPS = 0).

NDV: The number of times the velocity-space diagnostics (velocity distributions, spatial profiles, and phase space plots) are per- formed before they are displayed. In other words, these diag- nostics are displayed every NTV*NDV electron time steps. If NDV = 0, they are never displayed, but the velocity distribu- tions and spatial profiles are still calculated and stored to disk.

NPLOT: The number of plots which appear on the terminal at one time. This is useful for displaying multiple diagnostics simultane- ously. A value of NPLOT = 4 is reasonable.

MODED: The number of Fourier modes calculated by the ion density diagnostic.

NDD: The number of times the ion density diagnostic is performed be- fore the ion density in real space is displayed. In other words, the ion density is displayed every NTD*NDD electron time steps. If NDD = 0, the ion density is never displayed, but is still calculated and stored to disk.

The meanings of the variables in the third (extras) subgrouping are as follows:

where N0 is determined by the condition that NPX particles are are distributed in the region between EDGE < x < NX - EDGE, K0 = 2*PI*AMODEN/NX, and PHASE = FREQN*T0.

ANSE: The weight given the sinusoidal part of the preceding initial electron distribution.

where N0 is determined by the condition that NPX particles are are distributed in the region between EDGE < x < NX - EDGE, K0 = 2*PI*AMODEN/NX, and PHASE = FREQN*T0.

FREQN: The frequency which determines PHASE in the sinusoidal part of the preceding initial density distribution.

where N0 is determined by the condition that NPX particles are are distributed in the region between EDGE < x < NX - EDGE, K0 = 2*PI*AMODEN/NX, and PHASE = FREQN*T0.

ANSI: The weight given the sinusoidal part of the preceding initial ion distribution.

where K0 = 2*PI*AMODEX/NX.

FREQ: The frequency of the external pump field in the expression for for Ext(x) given above.

TRMP: The ramp-up time for the amplitude E1 and E2 of the external pump field in the expression for Ext(x) given above. The amp- litude E1 is determined by:

and the amplitude E2 is determined by:

E2 = BVR*(TIME/TRMP), if TIME < TRMP,
E2 = BVR, if TRMP < TIME < TOFF.

TOFF: The shut-off time for the amplitude E1 and E2 of the external pump field in the expression for Ext(x) given above. The amp- tude E1 = E2 = 0, if TIME > TOFF.

VTEST: Velocity of test particle, which is not allowed to change.

X0: Initial position of test particle.

1. INTRODUCTION

2. MAIN PART OF POST-PROCESSOR

The frequency spectrum P(w) of a complex function F(t) is defined to be:

where the sum is over the NT time points selected, and the user can choose the frequency values w over some range. The spectrum of the raw data is always calculated and displayed. As an option, the time history of the real and imaginary parts of each elected Fourier mode F(k,t), the evolution of the real and imaginary parts in the complex plane, and the time history of the logarithm of the modulus of each Fourier mode can be displayed.

Another option is to calculate the autocorrelation function C(tau) of each Fourier mode F(t), defined to be:

where G(t) = F(t) - F0, F0 = SUM?F(t)?/NT, and the sums are over the NT values of t. The frequency spectrum of the C(tau) will always be dis- played if this option is chosen. Furthermore, if the user had also cho- sen to display the time histories of the Fourier modes, then the time histories of the autocorrelation function will be similarly displayed.

1. INTRODUCTION

2. MAIN PART OF POST-PROCESSOR

where S(x) is the particle shape function, and the sum is over parti- cles. Of these, only the NPRO subset which was selected when the main simulation was run, was stored to disk and can be displayed. If at least the first three of these was chosen (i.e., if NPRS = 1 and NPRO was greater than or equal to 3), then it is also possible to display the following fluid quantities:

One can show that the fluid quantities are related to the velocity mo- ments according to the relations:

A time-average of the profiles can be optionally calculated, where

and the sum is over the NT time points selected, and x is the coordinate.

Finally, the scales of the plots can either be fixed (i.e., the same for all plots), or variable (i.e., allowing each plot to find its own best scale).

3. INPUT VARIABLES

Ti = T0 + DTC*LTS, where DTC = NTV*DT.

ITS: Increment between data points in the file to be used. This al- lows one to use only every other data point, for example. The time between data points is then DTI = DTC*ITS.

NTS: Total number of data points in the file to be used. The time corresponding to the last data point is Tf = Ti + DTI*(NTS - 1).

MTS: Number of data points used in time-average. The length of time the data is averaged over is then DTI*MTS. (MTS = 0 means that the time-average is not calculated.)

NDS: The number of spatial profiles to be displayed. NDS must be less than or equal to NPRO.

ION: Parameter which determines whether ion (ION = 1) or electron (ION = 0) quantities will be displayed. ION can be set to 1 only if MOVION was equal to 1.

IFL: Parameter which determines fluid variables (IFL = 1) or velocity moments (IFL = 0) will be displayed. IFL can be set to 1 only if NPRS = 1 and NPRO = 3 or 4.

NPLOT: The number of plots per page.

From these, nine executables and nine libraries must be created. They are:

There are, in addition, six alternative files:

from which one can create three executables and three libraries:

The command file contained in i.sh installs the complete system. To use it, type:

Eighteen jobs will be submitted to create the nine executables and nine libraries. Specific details regarding each of these files are in the sections following.

by executing the command file bldparm.sh This executable reads the simulation input file, input, allows the user to interactively update it, and then writes the updated file. More importantly, it also writes an include file, includa, and a namelist file, nlistf. The include file contains PARAMETER statements needed to dimension arrays in beps1.f, and the namelist file contains the remaining simulation input. These two files are the only way that data is passed to the simulation program.

Next, one must create seven graphics libraries. These are:

created by executing bldlibn.sh,

created by executing bldlibs.sh (which requires SunCore graphics),

created by executing bldlibt.sh,

created by executing nbldlibt.sh,

created by executing bldlibp.sh,

created by executing bldlibm.sh,

created by executing bldlibc.sh (which requires Ncargraphics).

In addition, one must create the simulation subroutine library

by executing bldsiml.sh.

Finally, one must create five executables for plotting or viewing graphics metafiles. These are:

created by executing bldgrms.sh (which requires SunCore graphics),

created by executing bldgrmt.sh,

created by executing nbldgrmt.sh,

created by executing bldgrmp.sh,

created by executing bldgrmc.sh (which requires Ncargraphics).

Alternative drivers for the Tektronix 4105 terminals exist using the commercial graphics packages DISSPLA and DI-3000. Specifically, if DISSPLA is available, the library nlibt.o can be replaced by libd.o, which is created by executing bldlibd.sh, and the executable ngramtk.out can be replaced by gramgt.out, which is created by executing bldgrmd.sh. Alternatively, if DI-3000 Graphics is available, then the graphics lib- rary nlibt.o can be replaced by libr.o, which is created by executing bldlibr.sh, and the executable ngramtk.out can be replaced by gramgr.out, which is created by executing bldgrmr.sh.

When all this is done, one executes the simulation program by typing:

This command file will compile beps1.f, using the include file includa. It will ask the user what graphics device to send the output to, and then will link with the appropriate graphics library as well as the simulation library simulib.o. If a metafile was selected, a dataset will be allocated for it, and defaults to agraph.m. The command file will then allocate the output datasets used by those diagnostics which the user had selected in the simulation input. The default names for the datasets are:

The simulation will then be executed, reading the namelist file for fur- ther input. If a metafile had been selected as the output device, the user is asked at the end of the simulation if the metafile should be plotted or viewed. If plotting was selected, the graphics is plotted by the command file contained in m.sh. If viewing was chosen, the graphics can be displayed on Tektronix or ascii terminals by executing the com- mand file contained in l.sh. It is also possible to plot metafiles independently of the main simulation program by typing:

where filename is the name of the metafile to be plotted. Similarly, the metafiles can be viewed independently of the main simulation by typing:

which is created by executing the command file bldanls.sh.

1. POTENTIAL DIAGNOSTIC

by executing the file bldcorp.sh. This executable reads the header of a file containing the potential data, and writes an include file, includa, which contains PARAMETER statements needed to dimension arrays in corre.f, according to the data in the header.

To execute this post-processor independently of the main simulation program, type:

where filename is the name of the file containing the potential data. This command file will compile corre.f using the include file includa, and link with the appropriate graphics library, as well as the post- processor library analysis.o. If a graphics metafile was selected, a dataset is allocated to contain it, which defaults to pgraph.m, and the user will be asked if it should be plotted or viewed.

2. DEBYE CLOUD DIAGNOSTIC

by executing the file bldcldp.sh. This executable reads the header of a file containing the potential data, and writes an include file, includa, which contains PARAMETER statements needed to dimension arrays in dcloud.f, according to the data in the header.

To execute this post-processor independently of the main simulation program, type:

where filesname is the name of the file containing the potential data and filername is the name of the file containing the potential data without the test charge (the reference run). This command file will compile dcloud.f using the include file includa, and link with the appropriate graphics library, as well as the post-processor library analysis.o. If a graphics metafile was selected, a dataset is alloca- ted to contain it, which defaults to cgraph.m, and the user will be asked if it should be plotted or viewed.

3. SPATIAL PROFILE DIAGNOSTIC

by executing the file bldpflp.sh. This executable reads the header of a file containing the profile data, and writes an include file, includa, which contains PARAMETER statements needed to dimension arrays in profile.f, according to the data in the header.

To execute this post-processor independently of the main simulation program, type:

where filename is the name of the file containing the profile data. This command file will compile profile.f using the include file includa, and link with the appropriate graphics library, as well as the post- processor library analysis.o. If a graphics metafile was selected, a dataset is allocated to contain it, which defaults to fgraph.m, and the user will be asked if it should be plotted or viewed.

These are just a guide, and you should explore at least one other set of parameters on your own. Note whether the agreement between simu- lation and theory gets better if the parameters you chose in the si- mulation more accurately reflect the approximations made in class.

B. Observe the results of the simulation. Here are few suggestions to to get you started. Feel free to watch other diagnostics as well.

a. Observe the evolution of phase space as the instability grows. Explain as many features as you can in terms of what you know from the Vlasov equation. Note the features that you cannot ex- plain or do not understand.

b. Describe the evolution of the potential. Can you correlate what happens to the potential plot with what is happening in phase space?

c. Examine the plots of total electric field, total kinetic energy, and net energy. Compare exchange of energy. How well is energy conserved?

d. Examine the plots of mode energy versus mode number. What is the fastest growing mode?

C. Run the potential post-processor. Set NTC=0 to turn off the auto- correlation calculation, since the instablity is not a stationary process in time. Can you estimate the frequency, growth rate, and wavenumber of the fastest growing instability (in physical units, i.e., in terms of plasma frequency, debye length, etc.)?

D. Write up a lab report which adequately addresses these questions. Include supporting evidence (plots, etc.). Due Feb. 24.

These are just a guide, and you could explore a different set of parameters, e.g., by using a different number of particles or system length.

B. Run the potential post-processor with the autocorrelation calcula- lation turned on (by using the default value of NTC). Observe both the correlated and uncorrelated potential as functions of time. What happens to the oscillations as the mode number increases? Can you estimate the frequency and damping rate? How well do they agree with the predictions for the Bohm-Gross waves?

C. Write up a brief lab report, and include supporting evidence. Due Mar. 11.

Field Energy History for mode 12

where w0**2 = wpe**2 + 3*(k0*vth)**2, and w0/k0 = 5*vth. Calculate the work done by the source and compare this with the expression (Ew**2/4*pi)*(Vph/(Vph - Vgr)), where Ew**2/8*pi is the wave energy, Vph is the phase velocity, and Vgr is the group velocity of the pump wave. Calculate the momentum given to the plasma and compare it with the expression (Ew**2/4*pi)/(Vph - Vgr). The following set of para- meters should get you started:

B. Repeat part A, where w0 and k0 are the same, but have the plasma moving with velocity v0 = 2*Vph. Verify that Ew**2/4*pi is the same as in part A. Show that the beam slows down and loses energy. Cal- culate the wave energy and momentum as in part A, and show that they satify the same relations as in part A. You can use the same para- meters as before, except let NPX = 0 and NPXB = 10240.

where w0**2 = wpe**2 + 3*(k0*vth)**2, and w0/k0 = 5*vth, but turn off the pump after 4 wave periods. Continue running the simulation for about 32 wave periods. Look for the parametric decay of this wave into back-scattered plasma waves and forward-going ion acoustic waves. Detect this by obtaining the power spectrum of the potential and ion density, using the appropriate post-processor, and distin- guishing between positive frequency (blue curve) and negative fre- quency (read curve) waves. (Hint: find the forward scattered ion wave first, then find the corresponding backscattered plasma wave.) Does k0 = ki + kr and w0 = wr + wi, where r and i refer to back- scattered and forward scattered waves, respectively? How does it agree with the dispersion relation? Find the growth rate of the ion acoustic wave. Below are some sample parameters to get you started:

B. Consider the same problem as in part A, but with fixed ions and look for Stimulated Compton Scattering by look at the power spectrum of the potential with the potential post-processor. Calculate the velo- city of the resonant particles from:

Remember kr is negative, check that this is true from your results. You may have to use the dispersion relation to find wr. Can you see any resonant particles in the phase space plots? You can use the same parameters as before, except that MOVION = 0.

B. Consider the same problem as in part A, but with k0*v0 = 2*Wpe and v0 = 4*vth.

COPYRIGHT 1988, REGENTS OF THE UNIVERSITY OF CALIFORNIA

Update: March 13, 1990