ngd.cpp

pic 模拟程序！面向对象

//--------------------------------------------------------------------
//
// File:	NGD.cpp
//
// Purpose:	implementations of the NGD (Neutral Gas Density) class
//
// Version:	$Id: ngd.cpp,v 1.30 2004/10/08 22:32:36 yew Exp $
//
// Copyright 2001, Tech-X Corporation, all rights reserved
//
//--------------------------------------------------------------------

#include "ngd.h"
#include "grid.h"

#ifdef _MSC_VER
#include <fstream>
using std::cout;
using std::cerr;
using std::endl;
#endif

const Scalar NGD::firstRevisionNumberWithNGD = 2.61;

NGD::NGD(Grid* pgrid, const ostring& gType, 
         const ostring& _analyticF, const Scalar& _gasDensity,  
         const Vector2& _p1Grid, const Vector2& _p2Grid, 
         const int& _discardDumpFileNGDDataFlag) 
  throw(Oops){

// Must create all memory or set pointer to zero at beginning

// Beginning of construction
  ptrLocalGrid = pgrid;
  gasType = gType;

  J = ptrLocalGrid->getJ(); 
  K = ptrLocalGrid->getK();
  NeutralGasDensity = new Scalar* [J];
  int j;
  int k;
  for (j = 0; j<J; j++) { 
    NeutralGasDensity[j] = new Scalar [K]; 
    for(k=0; k<K; k++) NeutralGasDensity[j][k]=0.0;
  }
  maxNGD = 0.0;

  ostring analyticF;
  // the following parser related lines are for the cross-platform compatability
#ifndef UNIX
  analyticF=_analyticF; 
  ostring endLine = "\n";
  analyticF=_analyticF + endLine;
#else
  analyticF=_analyticF + '\n';
#endif

// now set the Neutral Gas Density

  if( analyticF == (ostring)"0.0\n" ) { 

// this is the default value of analyticF
// for the case of uniform gasDensity
// calculated above and used here to
// initialize the NutralGasDensity
    try{
      setNeutralGasDensity( _gasDensity, _p1Grid, _p2Grid);
    }
    catch(Oops& oops){
      oops.prepend("NGD::NGD: Error: \n"); //SpatialRegion::initNGD
      throw oops;
    }

    maxNGD = _gasDensity;

  } else {
    /**
     * parse the function specified by the string contained in
     * analyticF to initialize the NeutralGasDensity and the
     * max value of the NeutralGasDensity.
     */
    try{
      setNeutralGasDensity(analyticF, _p1Grid, _p2Grid);
    }
    catch(Oops& oops){
      oops.prepend("NGD::NGD: Error: \n"); //OK
      throw oops;
    }
    
  }

  /**
   * set the boolean data member discardDumpFileNGDDataFlag 
   * according to the value passed from the input file
   */
  if ( !_discardDumpFileNGDDataFlag )
    discardDumpFileNGDDataFlag = false;
  else 
    discardDumpFileNGDDataFlag = true;
  
  /**
   * set TI_np2c to 0 as its default value
   */
  TI_np2c = 0;

  /**
   * set the excessNumIons pointer to 0 as its
   * default value. 
   */
  excessNumIons = 0;

  isTIOn = 0;

}

NGD::~NGD() {

  int j;
  for (j=0; j<J; j++) delete[] NeutralGasDensity[j];
  delete[] NeutralGasDensity;
  
  /**
   * dealocate if memory for np2c was allocated
   */
  if( TI_np2c ) {
    int j;
    for (j=0; j<J; j++) delete[] TI_np2c[j];
    delete[] TI_np2c;
  }  

// deallocate if memory for excessNumIons was allocated
  if( excessNumIons ) {
    int j;
    for (j=0; j<=J; j++) delete[] excessNumIons[j];
    delete[] excessNumIons;
  }   

}

/** 
 * helper function to check for out of bounds condition in both
 * the j and k indeces. For now the function will print which index
 * is out of bounds and in which function (given by char* place)
 * and exit to the system with status of 1. If the j, k are 
 * inbounds it will return true. dad 04/11/01.
 */
bool NGD::check_jk_bounds(const int& j, const int& k, char* place) 
  const throw(Oops){

  stringstream ss (stringstream::in | stringstream::out);

  if ( j < 0 || j > (J-1) ) {

    ss << "NGD::check_jk_bounds: Error: \n" << 
      "j = " << j << 
	    " is not in [0, " << J-1 << "]." << endl;

    std::string msg;
    ss >> msg;
    Oops oops(msg);
    throw oops;    // exit() NGD::getNGD:
   /*  NGD::getNGD_TI_np2c: 
       NGD::setNGD_TI_np2c: 
       NGD::setNGD_excessNumIons: 
       NGD::getNGD_excessNumIons:  
       NGD::add_to_NGD_excessNumIons:  
      */

  } else if ( k < 0 || k > (K-1) ) {

    ss << "NGD::check_jk_bounds: Error: \n" << 
      "k = " << k << 
	    " is not in [0, " << K-1 << "]." << endl;

    std::string msg;
    ss >> msg;
    Oops oops(msg);
    throw oops;    // exit()  NGD::getNGD:  
   /*   NGD::getNGD_TI_np2c:  
        NGD::setNGD_TI_np2c:  
        NGD::setNGD_excessNumIons:  
        NGD::getNGD_excessNumIons:  
        NGD::add_to_NGD_excessNumIons:  
      */

  }
  return true;
}

Scalar NGD::getNGD(const int& j, const int& k) const throw(Oops) {
  char* place = "NGD::getNGD(const int& j, const int& k)";

  try{
    check_jk_bounds(j, k, place); 
  }
  catch(Oops& oops){
    oops.prepend("NGD::getNGD: Error: \n");// lots of places
    throw oops;
  }

  return NeutralGasDensity[j][k];
}

//
// initializes the Neutral Gas Density
//
void NGD::setNeutralGasDensity(const Scalar& uniformGasDensity, 
                               const Vector2& p1Grid, const Vector2& p2Grid) 
  throw(Oops){
  int j,k;
  Vector2 x;
  Vector2 x1MKS = ptrLocalGrid->getMKS(p1Grid);
  Vector2 x2MKS = ptrLocalGrid->getMKS(p2Grid);
  
  

  if ( uniformGasDensity < 0.0 ) {
    stringstream ss (stringstream::in | stringstream::out);
    ss<< "NGD::setNeutralGasDensity: Error: \n"<<
      "Negative value for the uniform Neutral Gas Density" << 
	    endl << "This should never happen. Check the input file."<<endl;
  
    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops; // exit() "NGD::NGD: 
  }
  
  for (j=0; j < J; j++) {
    for (k=0; k < K; k++) { 
      x = ptrLocalGrid->getMKS(j,k);

      if ( x.e1() < x1MKS.e1() || x.e1() > x2MKS.e1() ||
	   x.e2() < x1MKS.e2() || x.e2() > x2MKS.e2() ) {
	NeutralGasDensity[j][k] = 0.0;
      } else {
	NeutralGasDensity[j][k] = uniformGasDensity;
      }
    }
  } 
}

void NGD::setNeutralGasDensity(const ostring & analyticF, 
                               const Vector2& p1Grid, const Vector2& p2Grid) 
  throw(Oops){
  int j,k;
  Vector2 x;
  Vector2 x1MKS = ptrLocalGrid->getMKS(p1Grid);
  Vector2 x2MKS = ptrLocalGrid->getMKS(p2Grid);

  stringstream ss (stringstream::in | stringstream::out);
  // maxNGD = 0.0; initial value set in the constructor
  
  for (j=0; j < J; j++) {
    for (k=0; k < K; k++){ 
      x = ptrLocalGrid->getMKS(j,k);
      if ( x.e1() < x1MKS.e1() || x.e1() > x2MKS.e1() ||
	   x.e2() < x1MKS.e2() || x.e2() > x2MKS.e2() )
	NeutralGasDensity[j][k] = 0.0;
      else {
	adv_eval->add_variable("x1",x.e1());
	adv_eval->add_variable("x2",x.e2());
	NeutralGasDensity[j][k] = adv_eval->Evaluate(analyticF.c_str());
	if ( NeutralGasDensity[j][k] < 0.0 ) {
    ss << "NGD::setNeutralGasDensity: Error: \n"<<
	    "Negative value for the Neutral Gas Density " <<
      endl << "while parsing the input function. This should never happen.";
	      
    
    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops;    // exit()  NGD::NGD
	  
	} else if ( NeutralGasDensity[j][k] > maxNGD ) {
	  maxNGD = NeutralGasDensity[j][k];
	}
      }
    }
  }
}

/**
 * get the number of physical particles per computational 
 * particle for a cell on the 2D grid. The cell is specified
 * by the j, k pair of indeces.
 */
Scalar NGD::getNGD_TI_np2c(const int& j, const int& k) const 
  throw(Oops){
  char* place = "NGD::getNGD_TI_np2c(const int& j, const int& k) const";
  try{
    check_jk_bounds(j, k, place);
  }
  catch(Oops& oops){
    oops.prepend("NGD::getNGD_TI_np2c: Error: \n");//OK
    throw oops;
  }

  Scalar result;
  stringstream ss (stringstream::in | stringstream::out);

  if ( TI_np2c ) 
    result = TI_np2c[j][k];
  else {
    ss<< "NGD::getNGD_TI_np2c: Error: \n" << 
      "Trying to access TI_np2c[" << j << "][" << k << "]\n" <<
      "while TI_np2c is set to 0. " << endl;
    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops;    // exit() CsMCTI::tunnelIonizeCs: 
      /*HeMCTI::tunnelIonizeHe:
      LiMCTI::tunnelIonize: 
      HMCTI::tunnelIonize: 
      MapNGDs::sendNGDsAndShift: 
      SpatialRegionBoundary::sendNGDStripe: 
      */
  }
  return result;
}

/**
 * set the number of physical particles per computational 
 * particle for a cell on the 2D grid. The cell is specified
 * by the j, k pair of indices.
 */
void NGD::setNGD_TI_np2c(const int j, const int k, Scalar value) 
  throw(Oops){
  stringstream ss (stringstream::in | stringstream::out);
  char* place = "NGD::setNGD_TI_np2c(const int j, const int k, Scalar value)";
  try{
    check_jk_bounds(j, k, place);  
  }
  catch(Oops& oops){
    oops.prepend("NGD::setNGD_TI_np2c: Error: \n"); //OK
    throw oops;
  }

  if ( TI_np2c ) 
    TI_np2c[j][k] = value;
  else {
    ss <<  "NGD::setNGD_TI_np2c: Error: \n" <<
        "Trying to access TI_np2c[" << j << "][" << k << "]"
         << endl << "while TI_np2c is set to 0. " << endl;
    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops;    // exit() MapNGDs::sendNGDsAndShift:  
    /*SpatialRegionBoundary::sendNGDStripe:
      SpatialRegionBoundary::recvNGDStripe:*/
  } 
}

/**
 * allocate memory for the TI_np2c and initialize it.
 */
void NGD::initNGD_TI_np2c(const Scalar& TI_numMacroParticlesPerCell, 
                          Grid* grid) 
  throw(Oops){
  if ( !TI_np2c ) {

    int j, k;
    Scalar cellVol;
    Scalar numCellNeutralAtoms;

    TI_np2c = new Scalar* [J];
    for (j = 0; j < J; j++) {
      TI_np2c[j] = new Scalar [K]; 
      for (k=0; k < K; k++) { 
        cellVol =  grid->cellVolume(j, k);
        try{
          numCellNeutralAtoms =  getNGD(j, k) * cellVol;
        }
        catch(Oops& oops2){
          oops2.prepend("NGD::initNGD_TI_np2c: Error: \n");//done
          throw oops2;
        }
        TI_np2c[j][k] = numCellNeutralAtoms/TI_numMacroParticlesPerCell;
      }
    }
  }
}

/**
 * allocate memory for excessNumIons and set its elements to the
 * default values of 0.0
 */

void NGD::initNGD_excessNumIons() {
  if ( !excessNumIons ) {
    int j, k;
    
    excessNumIons = new Scalar* [J+1];
    for (j = 0; j<=J; j++) {
      excessNumIons[j] = new Scalar [K+1];
      for ( k = 0; k <= K; k++ ) 
        excessNumIons[j][k] = 0.0;
    }
  }
}

/**
 * set the excess number of ions for cell (j, k)
 */
void NGD::setNGD_excessNumIons( const int& j, const int& k, 
                                const Scalar& value) 
  throw(Oops){
  stringstream ss (stringstream::in | stringstream::out);
  char* place = 
    "NGD::setNGD_excessNumIons(const int&, const int&, const Scalar&)";
  try{
    check_jk_bounds(j, k, place); 
  }
  catch(Oops& oops){
    oops.prepend("NGD::setNGD_excessNumIons: Error: \n"); //OK
    throw oops;
  }

  
  if ( excessNumIons ) 
    excessNumIons[j][k] = value;
  else {
    ss << "NGD::setNGD_excessNumIons: Error: \n"<<
      "Trying to access excessNumIons[" << j << "][" << k << "]" 
      << endl << "while excessNumIons is set to 0." << endl;

    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops;    // exit() CsMCTI::tunnelIonizeCs: Error:
    /*MapNGDs::sendNGDsAndShift
HeMCTI::tunnelIonizeHe: Error:
LiMCTI::tunnelIonize: Error: 
HMCTI::tunnelIonize: Error: 
SpatialRegionBoundary::recvNGDStripe*/

  }
}


/**
 * get the excess number of ions for cell (j, k)
 */
Scalar NGD::getNGD_excessNumIons( const int& j, const int& k) 
  const throw(Oops){
  stringstream ss (stringstream::in | stringstream::out);
  char* place = 
    "NGD::getNGD_excessNumIons( const int& j, const int& k)";
  try{
    check_jk_bounds(j, k, place);  
  }
  catch(Oops& oops){
    oops.prepend("NGD::getNGD_excessNumIons: Error: \n"); //OK
    throw oops;
  }

  Scalar result;
  if ( excessNumIons ) 
    result = excessNumIons[j][k];
  else {
    ss   <<"NGD::getNGD_excessNumIons: Error: \n" 
      "Trying to access excessNumIons[" << j << "][" << k << "] " 
      << endl << "while excessNumIons is set to 0." << endl;
    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops;    // exit()CsMCTI::tunnelIonizeCs: Error:
    /*MapNGDs::sendNGDsAndShift
HeMCTI::tunnelIonizeHe: Error:
LiMCTI::tunnelIonize: Error: 
HMCTI::tunnelIonize: Error: 
SpatialRegionBoundary::sendNGDStripe*/
  }
  return result;
}

/** 
 * add to the excess number of ions for cell (j, k)
 */
void NGD::add_to_NGD_excessNumIons( const int& j, const int& k, 
                                    const Scalar& value) 
  throw(Oops){
  stringstream ss (stringstream::in | stringstream::out);
  char* place = 
    "NGD::add_to_NGD_excessNumIons(const int&, const int, const Scalar&)";
  try{
    check_jk_bounds(j, k, place); 
  }
  catch(Oops& oops){
    oops.prepend("NGD::add_to_NGD_excessNumIons: Error: \n"); //OK
    throw oops;
  }
  if ( excessNumIons ) 
    excessNumIons[j][k] += value; 
  else {
    ss << "NGD::add_to_NGD_excessNumIons: Errors: \n"<<
      "Trying to access excessNumIons[" << j << "][" << k << "] " 
      << endl << "while excessNumIons is set to 0." << endl;
    std::string msg;
    ss>>msg;
    Oops oops(msg);
    throw oops;    // exit()// exit()CsMCTI::tunnelIonizeCs: Error:
    /*HeMCTI::tunnelIonizeHe: Error:
    LiMCTI::tunnelIonize: Error: 
    HMCTI::tunnelIonize: Error: */
  }
   
}

/**
 * Dump the data structures for the neutral gas density
 */

#ifdef HAVE_HDF5
int NGD::dumpH5(dumpHDF5 &dumpObj,int number)
{
  Vector2 **X = ptrLocalGrid->getX();
//  Scalar x,y;
  //cerr << "entered NGD::dumpH5.\n";


  hid_t fileId, groupId, datasetId,dataspaceId; //subgroupId ;
  herr_t status;
  hsize_t rank; 
  hsize_t     dimsf[3];   
  hsize_t     dims;
//  int         attrData;
  hid_t  attrdataspaceId,attributeId;
  hid_t scalarType = dumpObj.getHDF5ScalarType();
  string datasetName;

  //  strstream s;
#ifdef HAVE_SSTREAM
  stringstream s;
#else
  strstream s;