csmcti.cpp

pic 模拟程序！面向对象

#include <fstream>
#include <cstdio>

#include "CsMCTI.h"
#include "species.h"
#include "particle.h"
#include "ptclgrp.h"
#include "misc.h"

using namespace std;

//
// Tunneling ionization of Cs for only two of its electrons,
// the outer most ones, is implemented and thus Ztot = 2.  
// The ionization is from the ground states. 
//  
int const CsMCTI::Ztot = 2; 
                            

CsMCTI::CsMCTI(Scalar p, Scalar temp, Species* eSpecies, Species* iSpecies, 
               Species* _iSpeciesPlusPlus,
               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& numMacroParticlesPerCell, 
               const int& ETIPolarizationFlag, Scalar argEfieldFrequency,
               const int& discardDumpFileNGDDataFlag) throw (Oops)
  : MCTIPackage(Cs, (ostring)"Cs", p, temp, eSpecies, iSpecies, sList, dt, 
                ionzFraction, ecxFactor, icxFactor, _SR, _Min1MKS, _Max1MKS, 
                _Min2MKS, _Max2MKS, _delayTime, _stopTime, analyticF, 
                numMacroParticlesPerCell, argEfieldFrequency, 
                discardDumpFileNGDDataFlag){

  cout << endl << endl << "Creating CsMCTI" << endl;
  iSpeciesPlusPlus = _iSpeciesPlusPlus;
  nStar = new Scalar[Ztot]; // effective principal quantum #
  Zti   = new int[Ztot];    // charge of the atomic residue after each ti
  Ei    = new Scalar[Ztot]; // electron bound level ground state energies in au
  E0    = new Scalar[Ztot]; // auxiliary parameters 

  //
  // a constant factor for each bound state's
  // tunneling ionization probability rate
  // for linearly polarized E field  
  // 
  wLinearFactor = new Scalar[Ztot]; 

  // 
  // a constant factor for each bound state's    
  // tunneling ionization probability rate       
  // for circularly polarized and static E fields                         
  // 
  wCircularFactor = new Scalar[Ztot];                         

  //
  // memory for the max and min values of the external E field with
  // respect to the applicability of the AKD theory for
  // tunneling ionization 
  // 
  EmaxTI = new Scalar[Ztot];
  EminTI = new Scalar[Ztot];
 
  cout << "ATOMIC_ENERGY_UNIT_EV = " << ATOMIC_ENERGY_UNIT_EV << endl;
  Ei[0] = 3.8939/ATOMIC_ENERGY_UNIT_EV;
  Ei[1] = 25.1/ATOMIC_ENERGY_UNIT_EV;

  //
  // set the function pointer to the proper probability function
  // depending on the polarization of the electric field
  // 
  if ( ETIPolarizationFlag == 0 ) { // linear polarization
    fp = &CsMCTI::getLinearFieldProbabilityTI;
  } else if ( ETIPolarizationFlag == 1 ) { // circular polarization
    fp = &CsMCTI::getCircularFieldProbabilityTI;
  } else { // unrecognized value for the ETIPolarizationFlag
    stringstream ss (stringstream::in | stringstream::out);

    
    ss << "CsMCTI::CsMCTI: Error: \n" << 
         "Unrecognized value ETIPolarizationFlag = " << ETIPolarizationFlag 
         << endl << "There are currently defined valid values of 0 and 1\n" <<
         "Check your input paramter for ETIPolarizationFlag."<<endl;
         
    std::string msg;
    ss >> msg;
    Oops oops(msg);
    throw oops;    // exit() MCCParams::CreateCounterPart

  }
  Scalar tmp_dt = region->get_dt()/ATOMIC_TIME_UNIT;

  for ( int i = 0; i < Ztot; i++ ) {
    Zti[i]   = i + 1;
    nStar[i] = get_nStar(Ei[i], Zti[i]);
    E0[i] = pow(2.0*Ei[i], 1.5);
    wCircularFactor[i]  = pow((Scalar)4.0, (Scalar)nStar[i])*Ei[i]/nStar[i];
    wCircularFactor[i] /= sqrt(2.*M_PI)*exp(-2.0*nStar[i])*pow((Scalar)(2.0*nStar[i]), (Scalar)(2.0*nStar[i]-0.5));
    wCircularFactor[i] *= pow((Scalar)(2.0*E0[i]), (Scalar)(2.0*nStar[i]-1.0));
    if ( i == 1 ) // this is needed for l = 1, m = -1 for Cs+ -> Cs++
      wCircularFactor[1] *= (3./(2.*E0[i]));
    wLinearFactor[i]    = wCircularFactor[i]*sqrt(3.0/(M_PI*E0[i]));

    //
    // set the max field as the value at the inflection point of the 
    // probability rate function. 
    // 
    if ( i == 0 ) {
      if ( (2.0*nStar[i]-1.) > 1. ) {
        // 
        // There is an extremum point, so set the Emax to it. 
        // The max field should really be determined from
        // solving the time dependent Schrodinger equation. 
        // This needs further work...!
        // 
        EmaxTI[i] = 2.0*E0[i]/(3.0*(2.0*nStar[i]-1.)); 
      } else {
        stringstream ss (stringstream::in | stringstream::out);
        ss << "CsMCTI::CsMCTI: Error: \n" << 
              "WARNING: The tunneling ionziation probability rate \n" <<
              "does not have an extremum for the ionization of Cs. \n" <<
              "The max field for which the tunneling ionization \n"   << 
              "formula is applicable and the ionization rate for \n"  << 
              "larger fields is to be determined more precisely...\n";
              
        std::string msg;
        ss >> msg;
        Oops oops(msg);
        throw oops;    // exit()CreatecounterPart.
      }
    } else if ( i == 1 ) {
      if ( (2.0*nStar[i]-1.) > 1. ) {
        // 
        // There is an extremum point, so set the Emax to it. 
        // The max field should really be determined from
        // solving the time dependent Schrodinger equation. 
        // This needs further work...!
        // 
        EmaxTI[i] = 2.0*E0[i]/(3.0*(2.0*nStar[i]-2.)); 
      } else {
        stringstream ss (stringstream::in | stringstream::out);
        ss << "CsMCTI::CsMCTI: Error: \n" << 
           "WARNING: The tunneling ionziation probability rate " << endl
             << "does not have an extremum for the ionization of Cs." << endl
             << "The max field for which the tunneling ionization "   << endl
             << "formula is applicable and the ionization rate for "  << endl
             << "larger fields is to be determined more precisely..." << endl;
             

        std::string msg;
   	    ss >> msg;
   	    Oops oops(msg);
   	    throw oops;    // exit() CreatecounterPart.

      }
    }
    //
    // Check if the probability at EmaxTI and the time step of the simulation
    // is greater than one. If yes, print a warning. 
    // 
    Scalar temp = 0.0;
    try {
      temp = (this->*fp)(EmaxTI[i], tmp_dt, Zti[i]-1);
    }
    catch(Oops& oops){
      oops.prepend("CsMCTI::CsMCTI: Error:\n"); // OK
      throw oops;
    }

    if ( temp < 1.0 ) 
      cout << endl << endl << "WARNING: The probability at the max E field: " << endl 
           << "P(Emax[" << i << "] = " << EmaxTI[i] << ", dt = " << tmp_dt << ") = " 
           << temp << endl 
           << " is less than 1. The time step of the simulations is probably" << endl
           << "out of the regime of applicability of the tunneling ionization"
           << "theory!?" << endl << endl;
    
    cout << "Ei[" << i << "] = " << Ei[i] << endl
         << "Zti[" << i << "] = " << Zti[i] << endl
         << "nStar[" << i << "] = " << nStar[i] << endl
         << "E0[" << i << "] = " << E0[i] << endl << endl;

    cout << "P(Emax[" << i << "] = " << EmaxTI[i] << ", dt = " << tmp_dt << ") = " 
         << temp << endl << endl;

  /**
   * Set the value of the E field (in atomic units) below which 
   * the Keldish condition is violated. The probability for 
   * tunneling ionization for fields smaller than this value
   * will be set to zero. 
   * This value is given by w*sqrt(2*Ei) where w is the frequency
   * of the field, Ei is the energy of the bound state, both
   * quantities are in atomic units. Note that for Hydgrogen in
   * the ground state 2*Ei = 1. Therefore: 
   */
  EminTI[i] = EfieldFrequency*ATOMIC_TIME_UNIT*sqrt(2.*Ei[i]);
  cout << "EminTI[" << i << "] = " << EminTI[i] << endl;
  cout << "P(EminTI[" << i << "] = " << EminTI[i] << ", dt = " << tmp_dt 
       << ") = " << temp << endl; 
  }
}

CsMCTI::~CsMCTI() {
  if (nStar)
    delete [] nStar;
  if (Zti)
    delete [] Zti;
  if (Ei)
    delete [] Ei;
  if (E0)
    delete [] E0;
  if (wLinearFactor)
    delete [] wLinearFactor;
  if (wCircularFactor) 
    delete [] wCircularFactor;
  if (EmaxTI)
    delete [] EmaxTI;
  if (EminTI)
    delete [] EminTI;
}
/**
 * This member function manages the tunneling ionization of the
 * Cs atoms and Cs+ ions. 
 */
void CsMCTI::tunnelIonize(ParticleGroupList** pgList, ParticleList& pList) 
  throw(Oops){
  // 
  // attempt to tunnel ionize first Cs atoms for the current values of the 
  // E field: 
  // 
    try{
      tunnelIonizeCs(pList);
    }
    catch(Oops& oops2){
      oops2.prepend("CsMCTI::tunnelIonize: Error: \n");//done
      throw oops2;
    }

  // 
  // find the ParticleGroupList for the CsPlus species (ions) and then
  // attempt to tunnel ionize the particles in this list: 
  // 
  SpeciesList* speciesList = region->getSpeciesList();
  oopicListIter<Species> siter(*speciesList);
  for (siter.restart(); !siter.Done(); siter++) {
    Species* species = siter.current();
    int i = species->getID();
    //
    // we kept a pointer to the iSpecies in the BaseMCPackage class
    // and it points to the singly ionized atoms of this type, i.e.
    // Cs. 
    //
    if ( i == iSpecies->getID() ) {
      // 
      // now attempt to tunnel ionize the CsPlus ions
      // 
      try{
        tunnelIonizeCsPlus(*pgList[i], pList); 
      }
      catch(Oops& oops3){
        oops3.prepend("CsMCTI::tunnelIonize: Error: \n"); //done
        throw oops3;
      }

      break;
    }
  }
}

void CsMCTI::tunnelIonizeCs(ParticleList& pList) 
  throw(Oops){
  /**
   * loop over all cells of the grid and check for nonzero gas density
   * on each cell. If the tunneling ionization eoes not run fast enough, 
   * one way to speed it up could be to keep a list of the cells that have 
   * nonzero gas density and loop only over them. 
   * 
   * The J,K that we loop over are extracted from the Grid and should be
   * the same J,K that are data members of both the BaseMCPackage and NGD 
   * classes. 
   */
  
  int j,k;
  Vector2 x;      // temporary vector variable for positions on the grid
  Vector3 u;      // temporary vector variables for velocities of 
                  // created particles. The default value is 0. The effect of 
                  // temperature on the gas, that may lead to nonzero initial 
                  // velocities of created ions is not considered yet. 
  Vector3 Ecc;    // vector for the electric field at a grid cell center
  Scalar E;       // holds the magnitude of the electric field in atomic units
  Scalar dt;      // time interval during which to calculate the probability 
                  // for tunneling ionization
  Scalar cellVol; // volume of a cell on a grid
  Scalar numCellNeutralAtoms; // number of neutral atoms in a cell 
  int numMacroParticles;      // temporary varible for the number of macro
                              // particles to be created
  Scalar numIons; // temporary variable for the number of ions to create
  Scalar TI_np2c; // temporary variable to hold the number of physical ions
                  // in a macro parcticle for a given cell 
  Scalar expectedNumIons;  // temporary variable for the expected number of 
                           // atoms to become ions
  Scalar excessNumIons;    // temporary variable for the excess number of
                           // ions of a given grid cell
  static bool IonizationEnergyWarningFlag = false;

  /** 
   * set the default value for the velocity of created ion 
   * and electron particles
   */
  u.set_e1( 0.0 );
  u.set_e2( 0.0 );