mcc.cpp

pic 模拟程序！面向对象

/*
=====================================================================
mcc.cpp

Methods for the MCC (Monte Carlo collision) class.

0.9	(JohnV 02-06-95) Initial draft.
1.01  (JohnV 01-02-97) Restructure for flexibility.
1.1   (JohnV 08-04-97) Added MCCPackages for Ar, Ne and Xe.
1.11  (JohnV 08-15-97) Repaired disparity between cross sections and
      thresholds for inelastic collisions.
1.12  (JohnV 10-01-97) Fixed mcc to work with supercycled timesteps.
1.13  (KC 10-20-97) Fixed mcc to work with M/m < 2.
1.14  (JohnV 02-11-99) Added He package (from Y. Ikeda)
2.5   (Bruhwiler 03-03-00) Added Lithium package (from S.M. Younger)
2.51  (Bruhwiler 08-29-00) Modified Lithium package to be relativistic
2.52  (Bruhwiler 09-03-00) added relativistic elastic scattering for
       Li -- this required changes in CrossSection and FuncCrossSection
2.53  (Bruhwiler 09-26-00) added capability to limit MCC region (i.e.
       limit the neutral gas) to a subregion of the grid.
2.54  (JohnV 08NOV2000) fixed glitch in e-Xe elastic scattering cross section.

Consider how to make cross sections work for adjusted-mass particles
		
=====================================================================
*/
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <cstdio>
#include <txc_streams.h>

#ifdef _MSC_VER
#include <fstream>
#else
//#include <txc_streams.h>
using namespace std;
#endif

#include "mcc.h"
#include "species.h"
#include "maxwelln.h"
#include "particle.h"
#include "ptclgrp.h"
#include "misc.h"

const Scalar	NeVperTORR =	8.3221e20;
const Scalar	MAX_ENERGY = 1000;

#define pow(x,y) pow((Scalar)(x),(Scalar)(y))

//-------------------------------------------------------------------
// CrossSection: provides basic cross section handling capability


//-------------------------------------------------------------------
//	eCrossSection methods

eCrossSection::eCrossSection(Scalar e0, Scalar e1, Scalar e2, Scalar sMax, Type _type)
{
	energy0 = e0;
	energy1 = e1;
	energy2 = e2;
	sigmaMax = sMax;
	const1 = sigmaMax/(energy1 - energy0);
	const3 = sigmaMax*energy2/log(energy2);
	type = _type;
}

Scalar eCrossSection::sigma(Scalar energy)
{
	if (energy <= energy0) return 0;
	if (energy <= energy1) return (energy - energy0)*const1;
	if (energy <= energy2) return sigmaMax;
	return const3*log(energy)/energy;
}


//-------------------------------------------------------------------
// Argon ionization cross section (Rapp & Golden) - energies in eV

Scalar argonIoniz::sigma(Scalar energy)
{										  // 
   if (energy < 15.76) return 0;
	else if (energy < 79) 
	  return 1.7155e-18*(energy-15.76)/(energy*energy)*log(0.0634*energy);
	else return 2.648e-18*(energy-15.76)/(energy*energy)*log(0.0344*energy);
}

//-------------------------------------------------------------------
//	iCrossSection methods

iCrossSection::iCrossSection(Scalar _a, Scalar _b, Type _type)
{
	a = _a;
	b = _b;
	type = _type;
}

Scalar iCrossSection::sigma(Scalar energy)
{
	return a + b/(sqrt(energy) + 1E-30);
}

//--------------------------------------------------------------------
// MCC: flexible collision package for handling multiple gases

MCC::~MCC()
{
  packageList.deleteAll(); // remove all list elements and data for packages
  pList.deleteAll();       // delete the list elements and data for the particles
                               
}

ParticleList& MCC::collide(ParticleGroupList& pgList)
  throw(Oops){
  oopicListIter<MCCPackage> packageIter(packageList);
  for (packageIter.restart(); !packageIter.Done(); packageIter++){
    try{
	    packageIter()->collide(pgList, pList);
    }
    catch(Oops& oops3){
      oops3.prepend("MCC::collide: Error: \n");
      throw oops3;
    }
  }
  return pList;
}

/**
 * MCTI: Monte Carlo Tunneling Ionization package methods 
 * for handling multiple gases
 * dad 03/27/01
 */
MCTI::~MCTI() {
  packageList.deleteAll(); // remove all list elements and data for the MCTI packages
  pList.deleteAll();       // delete the list elements and data for the particles
}

void MCTI::tunnelIonize(ParticleGroupList** pgList)
throw(Oops){
  oopicListIter<MCTIPackage> packageIter(packageList);
  for (packageIter.restart(); !packageIter.Done(); packageIter++){
    try{
	    packageIter()->tunnelIonize(pgList, pList);
    }
    catch(Oops& oops2){
      oops2.prepend("MCTI::tunnelIonize: Error: \n");//done
      throw oops2;
    }
  }
}

// BaseMCPackge
BaseMCPackage::BaseMCPackage(GAS_TYPE _type, const ostring& strGasType, 
                             const Scalar& temperature, 
                             Species* _eSpecies, Species* _iSpecies,
                             const Scalar& pressure, SpatialRegion* _SR,
                             Scalar _Min1MKS, Scalar _Max1MKS, 
                             Scalar _Min2MKS, Scalar _Max2MKS,
                             const ostring& _analyticF, 
                             const int& discardDumpFileNGDDataFlag)
  throw(Oops){
  type = _type;
  gasTypeString = strGasType;
  temp = temperature;
  gasDensity = NeVperTORR * pressure / ( temp + 1.e-30 );

  //
  // set default values for the following two pointers
  // (they are used as handles)
  //
  pList = 0;
  pg = 0;

  // initialize the Species pointers
  eSpecies = _eSpecies;
  iSpecies = _iSpecies;

  // gas may be confined to a portion of the grid (Bruhwiler 09/25/2000)
  Min1MKS = MIN(_Min1MKS,_Max1MKS);
  Max1MKS = MAX(_Max1MKS,_Min1MKS);
  Min2MKS = MIN(_Min2MKS,_Max2MKS);
  Max2MKS = MAX(_Max2MKS,_Min2MKS);

  region = _SR;
  grid   = region -> get_grid();
  J      = grid -> getJ();
  K      = grid -> getK();
  
  // skip the following logic in default case
  if (_Min1MKS > 0. && _Max1MKS > 0. && _Min2MKS > 0. && _Max2MKS > 0.) {
    // determine the maxima for this grid.  Clip the load if it's too large.
    Scalar x1min = grid->getX()[0][0].e1();
    Scalar x2min = grid->getX()[0][0].e2();
    Scalar x1max = grid->getX()[grid->getJ()][0].e1();
    Scalar x2max = grid->getX()[0][grid->getK()].e2();

    Min1MKS = MAX(x1min,Min1MKS);
    Max1MKS = MIN(x1max,Max1MKS);
    Min2MKS = MAX(x2min,Min2MKS);
    Max2MKS = MIN(x2max,Max2MKS);

    p1 = Vector2(Min1MKS,Min2MKS);
    p2 = Vector2(Max1MKS,Max2MKS);
    deltaP = p2-p1;
    
    p1Grid = grid->getGridCoords(p1);
    p2Grid = grid->getGridCoords(p2)+Vector2(1,1);
    Vector2 xt = grid->getGridCoords(p2);

    if ( J <= (int)p2Grid.e1() )  p2Grid.set_e1( (Scalar)J );
    if ( K <= (int)p2Grid.e2() )  p2Grid.set_e2( (Scalar)K );

    if (grid->query_geometry() == ZXGEOM) rz = FALSE;
    else if (grid->query_geometry() == ZRGEOM) rz = TRUE;
    else {
      stringstream ss (stringstream::in | stringstream::out);
      ss << "BaseMCPackage::BaseMCPackage: Error: \n"<<
        "MCCPackage constructor doesn't recognize the geometery flag" << TXSTD::endl;

      std::string msg;
      ss >> msg;
      Oops oops(msg);
      throw oops;    // exit() not called 

    }
  }  // end of logic regarding limited region for gas density
  else {
    p1Grid.set_e1( 0. );
    p1Grid.set_e2( 0. );
    p2Grid.set_e1( (Scalar)J );
    p2Grid.set_e2( (Scalar)K );
  }

  /**
   * at this point I am ready to initialize the NGD.
   * pass the string (ostring) of the function to be used for the evaluation of
   * the NeutralGasDensity data member of the Grid class. 
   * dad: Wed Nov 15 14:39:43 MST 2000}
   */
  try{
    ptrNGD = region->initNGD(gasTypeString, _analyticF, gasDensity, p1Grid, p2Grid, 
                           discardDumpFileNGDDataFlag);
  }
  catch(Oops& oops3){
    oops3.prepend("BaseMCPackage::BaseMCPackage: Error: \n");
    throw oops3;
  }
  maxNGD = ptrNGD->getMaxNGD();
  gasDensity = maxNGD;
  
}
//----------------------------------------------------------------
// MCCPackage: individual gas packages

MCCPackage::MCCPackage(GAS_TYPE _type, const ostring& strGasType, 
                       Scalar pressure, Scalar _temp,
                       Species* _eSpecies, Species* _iSpecies, SpeciesList& sList,
                       Scalar dt, int _ionzFraction, Scalar _ecxFactor, Scalar _icxFactor,
                       SpatialRegion* _SR, Scalar _Min1MKS, Scalar _Max1MKS,
                       Scalar _Min2MKS, Scalar _Max2MKS, Scalar _delayTime, Scalar _stopTime, 
                       const ostring& _analyticF, const int& discardDumpFileNGDDataFlag) 
  : BaseMCPackage( _type, strGasType, _temp, _eSpecies, _iSpecies, pressure, _SR, 
                   _Min1MKS, _Max1MKS, _Min2MKS, _Max2MKS, _analyticF, 
                   discardDumpFileNGDDataFlag)
{
  ionzFraction = _ionzFraction;
  ecxFactor = _ecxFactor;

  icxFactor = _icxFactor;
  neCX = niCX = 0;
  eCX = iCX = NULL;

  delayTime = _delayTime;
  stopTime  = _stopTime;

  addCrossSections();
}

void MCCPackage::init(SpeciesList& sList, Scalar dt)
{
  if (gasDensity > 0) {
    extra = new Scalar[sList.nItems()];
    collProb = new Scalar[sList.nItems()];
    // eEnergyScale = fabs(0.5/eSpecies->get_q_over_m());
    // iEnergyScale = fabs(0.5/iSpecies->get_q_over_m());
    eEnergyScale = fabs(0.5/PROTON_CHARGE*eSpecies->get_m());
    iEnergyScale = fabs(0.5/PROTON_CHARGE*iSpecies->get_m());
    
    // check for electron and ion collisions
    maxSigmaVe = maxSigmaVi = 0;
    oopicListIter<Species> sIter(sList);
    for (sIter.restart(); !sIter.Done(); sIter++){
      int id = sIter.current()->getID();
      extra[id] = 0.5;
      switch (sIter()->get_collisionModel()) {
      case none:
	break;
      case boltzmann:
	//			BoltzmannInit();
	break;
      case electron:
      case test:
	if (maxSigmaVe == 0) maxSigmaVe = sumSigmaVe();
	collProb[id] = 1 - exp(-gasDensity*maxSigmaVe*dt*sIter()->get_subcycleIndex()
			       /sIter()->get_supercycleIndex());
	break;
      case ion:
	if (maxSigmaVi == 0) maxSigmaVi = sumSigmaVi();
	collProb[id] = 1 - exp(-gasDensity*maxSigmaVi*dt*sIter()->get_subcycleIndex()
			       /sIter()->get_supercycleIndex());
	break;
      }
    }
    
    //	create isotropic maxwellian gas distribution with ion characteristics
    gasDistribution = new Maxwellian(iSpecies);
    gasDistribution->setThermal(temp, EV);
  }
  else extra = collProb = NULL;
}

MCCPackage::~MCCPackage()
{
  if (eCX) {
	 for (int i=0; i<neCX; i++) delete eCX[i];
	 delete[] eCX;
  }
  if (iCX) {
	 for (int i=0; i<niCX; i++) delete iCX[i];
	 delete[] iCX;
  }
  delete[] extra;
  delete[] collProb;
  delete gasDistribution;
}

void MCCPackage::addCrossSections()
{
  neCX = niCX = 0;
  eCX = iCX = NULL;
}

void MCCPackage::collide(ParticleGroupList& pgList, ParticleList& _pList)
  throw(Oops){
  if (!pgList.head || !collProb) return;
  
  // return if simulation time is not within the specified time interval
  Scalar time = region->getTime();
  //  cout << "_____________________" << endl;
  //  cout << " simulation time is: " << time << endl << endl;
  if ( time<delayTime ) return;
  if ( time>stopTime && stopTime>0.) return;
  
  Species* species = pgList.head->data->get_species();
  
  CollisionModel collisionModel = species->get_collisionModel();
  if (collisionModel == none) return;
  int id = species->getID();
  if (collProb[id] < 1e-20) return;
  
  pList = &_pList;
  int n = 0;
  oopicListIter<ParticleGroup> pgIter(pgList);
  for (pgIter.restart(); !pgIter.Done(); pgIter++) {
    // n is the total # of particles in all relevant particle groups
    n += pgIter()->get_n();
  }
  if (n==0) return; // return if no particles
  
  // number of potential collisions is n times the max. collision prob.

  ///cout << "n = " << n << endl;

  extra[id] += n * collProb[id];
  int nCollisions = (int)extra[id];
  extra[id] -= nCollisions;

  ///cout << "nCollisions = " << nCollisions << endl;
  
  // holds position vector of particle that might undergo collision