uspas01.inp

pic 模拟程序！面向对象

uspas01
{

This input file, uspas01.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 5 MeV electron bunch enters vacuum in 2-D Cartesian geometry --
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) How are the electrons injected into the simulation?
  b) What is the particle distribution in x, y, Ux, Uy?
  c) What sort of fields are generated by the electron bunch?
  d) How important are the boundary conditions?
  e) What happens if the time step violates the Courant condition?
}

// 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

// Next, define the parameters of the high-energy electron beam.
  beamEnergyEV = 5.0e+06

  totalNumBeam = 100.0e+6
  totalBeamCharge = totalNumBeam * electronCharge
  rmsBeamWidth  = 1.0e-04
  rmsBeamLength = 1.0e-04
  totalBeamWidth  = 6 * rmsBeamWidth 
  totalBeamLength = 6 * rmsBeamLength
  totalBeamArea = totalBeamWidth * totalBeamWidth
  rmsBeamVolume = rmsBeamWidth * rmsBeamWidth * rmsBeamLength

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

// Define the number of grids in X and Y
  numYgrids = 32
  numXgrids = 32
  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 = 10
  maxLengthMKS = lengthFactor * rmsBeamLength
  xGridSize = maxLengthMKS / numYgrids
  widthFactor = 10
  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.
}

// 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
//  rmsDiagnosticsFlag = 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
}

}