uspas03.inp

pic 模拟程序！面向对象

uspas03
{

This input file, uspas03.inp, was written by David Bruhwiler
for "Object Oriented Computational Accelerator Physics," a
two week course presented at the University of Colorado for
the U.S. Particle Accelerator School in June, 2001.

A 30 GeV electron bunch enters a plasma in 2-D Cartesian geometry --
The parameters are comparable to those in the E-157 proof-of-principle
  PWFA experiment at SLAC.
Low resolution and few particles are used here for testing purposes.

Moving window:
  Once the electron beam has entered the grid and is close to the far
  edge of the simulation region, a moving window algorithm is invoked
  so that the beam can be modeled for long times.

Boundary conditions:
  The simulation region is bounded by perfect conductors.  This is
  not optimal, but it correctly handles the issue of particles that
  exit the simulation.

Topics for discussion:
  a) Why do the plasma electrons get "blown out"?
  b) What is the cause of the large Ex behind the particle bunch?
  c) How could you use this system as a particle accelerator?
  d) What happens if the plasma density is increased?
}

// Define variables that can be used throughout this input file.

Variables
{
// First, define some useful constants.
  speedOfLight = 2.99792458e+08
  electronMass = 9.1093897e-31
  unitCharge = electronMass * 1.75881962e11
  electronCharge = -1. * unitCharge
  electronMassEV = electronMass * speedOfLight * speedOfLight / unitCharge
  ionCharge = unitCharge
  unitMassMKS = electronMass / 5.48579903e-04
  lithiumMassNum = 6.942
  lithiumMass = unitMassMKS * lithiumMassNum

// Next, define the parameters of the high-energy electron beam.
  beamEnergyEV = 30.0e+09
  beamGammaMin1 = beamEnergyEV / electronMassEV
  beamGamma = 1 + beamGammaMin1
  beamBetaGamma = sqrt( beamGammaMin1 * (beamGammaMin1+2) )
  beamBeta = beamBetaGamma / beamGamma

  totalNumBeam = 6.e+12
  totalBeamCharge = totalNumBeam * electronCharge
  rmsBeamWidth  = 1.0e-04
  rmsBeamLength = 5.0e-04
  rmsBeamTime = rmsBeamLength / speedOfLight
  totalBeamWidth  = 6 * rmsBeamWidth 
  totalBeamLength = 6 * rmsBeamLength
  totalBeamArea = totalBeamWidth * totalBeamWidth
  rmsBeamVolume = rmsBeamWidth * rmsBeamWidth * rmsBeamLength

  rmsNormalizedEmittance = 1.0e-07
  rmsThermalBeta = rmsNormalizedEmittance / rmsBeamWidth
  rmsThermalGamma = 1. / sqrt(1.-rmsThermalBeta*rmsThermalBeta)
  rmsVelocityEV = (rmsThermalGamma-1.)*electronMassEV

// Define the number of grids in X and Y
  numYgrids = 64
  numXgrids = 64
  numCells = numXgrids * numYgrids

// Calculate the size of the simulation region, grid spacings, time step.
// We are assuming the same grid size in both z and r	
  lengthFactor = 32
  maxLengthMKS = lengthFactor * rmsBeamLength
  xGridSize = maxLengthMKS / numYgrids
  widthFactor = 32
  maxWidthMKS = widthFactor * rmsBeamWidth
  yGridSize = maxWidthMKS / numYgrids
  effGridSize = 1. / sqrt( 1./(xGridSize*xGridSize) + 1./(yGridSize*yGridSize) )
  timeStep = 0.9 * effGridSize / speedOfLight

  yMiddle = 0.5 * maxWidthMKS
  numXgridsAcrossBeam = 6 * numXgrids / lengthFactor
  numYgridsAcrossBeam = 6 * numYgrids /  widthFactor

// Number of beam particles
  numBeamPtclsPerCell = 4
  numBeamCells = numXgridsAcrossBeam * numYgridsAcrossBeam
  numBeamPtcls = numBeamPtclsPerCell * numBeamCells
  beamNumRatio = totalNumBeam / numBeamPtcls

// Intermediate calculations for modeling Gaussian shape of the beam.
  invSigXsq = 1.0 / ( rmsBeamLength * rmsBeamLength )
  invSigYsq = 1.0 / ( rmsBeamWidth  * rmsBeamWidth  )
  invSigTsq = invSigXsq * speedOfLight * speedOfLight

// This is the desired delay time before the moving window algorithm activates.
  movingWindowDelay = 0.95 * maxLengthMKS / speedOfLight

// Calculate peak currents for defining emission of the high-energy beam.
  peakCurrentDensity=totalBeamCharge*speedOfLight/rmsBeamVolume
  peakCurrent = peakCurrentDensity * totalBeamArea
  pulseLengthSec = totalBeamLength / speedOfLight
  oneHalfPulse = pulseLengthSec/2.

// Define the plasma density, number of plasma electron macro-particles, etc.
  plasmaDensityMKS = 2.e20
  simulationVolume = maxWidthMKS * 1.0 * maxLengthMKS
  totalNumPlasma = plasmaDensityMKS * simulationVolume
  numPtclsPerCell = 4
  numPlasmaPtcls = numPtclsPerCell * numCells
  plasmaNumRatio = totalNumPlasma / numPlasmaPtcls
}

// This simulation has only one "region", which contains
//   the grid, all particles, etc.
Region
{

// Define the grid for this region.
Grid
{

// Define number of grids along X-axis and physical coordinates.
  J = numXgrids
  x1s = 0.0
  x1f = maxLengthMKS
  n1 = 1.0

// Define number of grids along Y-axis and physical coordinates.
  K = numYgrids
  x2s = 0.0
  x2f = maxWidthMKS
  n2 = 1.0

// Specify Cartesian geometry
  Geometry = 1
}

// Specify "control" parameters for this region
Control
{
// Specify the time step.
  dt = timeStep

// Turn on the moving window algorithm.
  movingWindow = 1
  shiftDelayTime = movingWindowDelay
}

// Define the beam electrons.
Species
{
  name = beam_electrons
  m = electronMass
  q = electronCharge
}

// Define the plasma ions.
Species
{
  name = plasma_ions
  m = lithiumMass
  q = ionCharge
}

// Load the plasma ions over the entire simulation region.
Load
{
  speciesName = plasma_ions
  density = plasmaDensityMKS
  x1MinMKS = 0.0
  x1MaxMKS = maxLengthMKS
  x2MinMKS = 0.0
  x2MaxMKS = maxWidthMKS

// This specifies a static uniform background (no macro-particles).
  np2c = 0
}

// Define the plasma electrons.
Species
{
  name = plasma_electrons
  m = electronMass
  q = electronCharge
}

// Load the plasma electrons over the entire simulation region, but
//   leave the last dz strip of cells empty, because this strip must
//   be handled separately to accomodate the moving window algorithm.
Load
{
  speciesName = plasma_electrons
  density = plasmaDensityMKS
  x1MinMKS = 0.0
  x1MaxMKS = maxLengthMKS - xGridSize
  x2MinMKS = 0.0
  x2MaxMKS = maxWidthMKS
  np2c = 2 * plasmaNumRatio

// Specify loading that is more uniform than random
  LoadMethodFlag = 1
}

// Load the plasma electrons into the last dz strip of cells, which was
//   omitted by the load instruction above.
Load
{
// Name this load group "shiftLoad" so that the moving window algorithm
//   knows to invoke it every time the simulation window is shifted.
  Name = shiftLoad
  speciesName = plasma_electrons
  density = plasmaDensityMKS
  x1MinMKS = maxLengthMKS - xGridSize
  x1MaxMKS = maxLengthMKS
  x2MinMKS = 0.0
  x2MaxMKS = maxWidthMKS
  np2c = 2 * plasmaNumRatio

// Specify loading that is more uniform than random
	LoadMethodFlag = 1
}

// Define the beam emitter, which introduces the high-energy beam into the
// simulation.
BeamEmitter
{
  speciesName = beam_electrons
  I = peakCurrent

// Define the 2-D function F(x,t) that specifies beam emission profile.
  xtFlag = 3
  nIntervals = 32
  F=exp(-invSigYsq*(x-yMiddle)*(x-yMiddle))*exp(-invSigTsq*(t-oneHalfPulse)*(t-oneHalfPulse))*step(pulseLengthSec-t)

// Macroparticles are emitted from the left boundary,
//   close to the axis of symmetry.
  j1 = 0
  j2 = 0
  k1 = (numYgrids - numYgridsAcrossBeam) / 2.
  k2 = (numYgrids + numYgridsAcrossBeam) / 2.
  normal = 1
  np2c = beamNumRatio

// Emit particles, directed along the Z-axis,  with specified
//    energy and temperature.
  units = EV
  v1drift = beamEnergyEV
  v1thermal = rmsVelocityEV
  v2thermal = rmsVelocityEV
  v3thermal = 0.
}

// Specify a perfect conductor along the left boundary.  This serves
//    as a particle boundary condition (catches particles that leave
//    the simulation) and as a field boundary condition (E_y is forced
//    to vanish).
Conductor
{
  j1 = 0
  j2 = 0
  k1 = 0
  k2 = numYgrids
  normal = 1
}

// Specify a perfect conductor along the top boundary.  This serves as a
//   particle boundary condition (catches particles that leave the simulation)
//   and as a field boundary condition (E_x is forced to vanish).
Conductor
{
  j1 = 0
  j2 = numXgrids
  k1 = numYgrids
  k2 = numYgrids
  normal = -1
}

// Specify a perfect conductor along the bottom boundary.  This serves as a
//   particle boundary condition (catches particles that leave the simulation)
//   and as a field boundary condition (E_x is forced to vanish).
Conductor
{
  j1 = 0
  j2 = numXgrids
  k1 = 0
  k2 = 0
  normal = 1
}

// Specify a perfect conductor along the right boundary.  This serves as a
//   particle boundary condition (catches particles that leave the simulation)
//   and as a field boundary condition (E_y is forced to vanish).
Conductor
{
  j1 = numXgrids
  j2 = numXgrids
  k1 = numYgrids
  k2 = 0
  normal = -1
}

}