csmcti.cpp

pic 模拟程序！面向对象

  u.set_e3( 0.0 );
   
  for ( j = 0; j < J; j++ ) {
    for ( k = 0; k < K; k++ ) {


      cellVol =  grid->cellVolume(j, k);
      try{
        numCellNeutralAtoms =  (ptrNGD->getNGD(j, k)) * cellVol;    
      }
      catch(Oops& oops2){
        oops2.prepend("CsMCTI::tunnelIonizeCs: Error: \n");//done
        throw oops2;
      }
      /**
       * proceed with the creation of macro particles of ions only
       * if the remaining number of neutral atoms is greater
       * than zero
       */
      if ( numCellNeutralAtoms < 1.0 )
        continue;
        
      /**
       * calculate the magnitude of the electric field at the 
       * cell's center in atomic units of field
       */
      x.set_e1( static_cast<Scalar>(j) + 0.5 );
      x.set_e2( static_cast<Scalar>(k) + 0.5 );
      Ecc = (region->get_fields())->E(x); 
      E   = Ecc.magnitude()/ATOMIC_E_FIELD_UNIT_MKS; 
      
      /**
       * calculate the time step in atomic units of time
       */
      dt  = region->get_dt() / ATOMIC_TIME_UNIT;
      
      /**
       * get the exptected number of ions to be created via 
       * the tunneling ionization. 
       *
       * Note that the probability function for ionization 
       * is coded specifically for tunneling from the ground state 
       * of Hydrogen. 
       */
      try{
        expectedNumIons = numCellNeutralAtoms * TIProbability(fp, E, dt, Zti[0]-1); 
      }
      catch(Oops& oops2){
        oops2.prepend("CsMCTI::tunnelIonizeCs: Error: \n"); //OK
        throw oops2;
      }
      //
      // is the energy of ionization bigger than 0.1% of the field energy
      // if yes print a warning but just once
      // 
      if ( !IonizationEnergyWarningFlag ) {
        if ( expectedNumIons > 1.0 ) 
          if ( expectedNumIons*0.5*ATOMIC_ENERGY_UNIT_EV/
               getCellEfieldEnergy(j, k) > 0.001 ) { 
            cerr << "WARNING: Detected Ionization Energy larger " << endl
                 << "than 0.1 % of the cell electric field energy!!!"
                 << endl;
            IonizationEnergyWarningFlag = true;
          }
      }
      /**
       * reduce the numCellNeutralAtoms by the expected number
       * of atoms to become ions and update the neutral gas
       * density for the cell
       */
      Scalar stmp = (numCellNeutralAtoms - expectedNumIons)/cellVol; 
      if ( stmp < 0.0 ) {
        cerr << "numCellNeutralAtoms - expectedNumIons = " 
             <<  numCellNeutralAtoms - expectedNumIons << endl
             << "is negative. This should never happen." << endl;
        stmp = 0.;
        expectedNumIons = numCellNeutralAtoms;
      }
      ptrNGD->set(j, k,  stmp);
      
      /** 
       * add to the excess number of ions for the (j, k)th cell
       */
      try{
        ptrNGD->add_to_NGD_excessNumIons(j, k, expectedNumIons);
        
        /**
        * create macro particles with ions only if 
        * the current's cell excessNumIons > TI_np2c
        */
        TI_np2c = ptrNGD->getNGD_TI_np2c(j,k);
        excessNumIons = ptrNGD->getNGD_excessNumIons(j, k);
      }
      catch(Oops& oops){
        oops.prepend("CsMCTI::tunnelIonizeCs: Error: \n"); //OK CsMCTI::tunnelIonize:
        throw oops;
      }


      /**
       * create the macro particles with ions only if 
       * the current excessNumIons is greater than
       * the number of physical ions per macro particle
       */
      if ( excessNumIons > TI_np2c ) {
        /** 
         * calculate the number of macro particles to create
         */
        numMacroParticles = static_cast<int>(excessNumIons/TI_np2c);
        /**
         * reset the excessNumIons[j][k] data member of class NGD
         */
        numIons = TI_np2c * (static_cast<Scalar>(numMacroParticles));
        excessNumIons -= numIons;
        try{
          ptrNGD->setNGD_excessNumIons(j, k, excessNumIons);
        }
        catch(Oops& oops){
          oops.prepend("CsMCTI::tunnelIonizeCs: Error: \n"); //OK
          throw oops;
        }
        /**
         * create the macro particles with 0 velocities and positions
         * picked randomly on the cell (in grid coordinates)
         */
        int MPcounter;
        for ( MPcounter = 0; MPcounter < numMacroParticles; MPcounter++ ) {
          /** 
           * create the particles for both ions and electrons
           * with random positions in the cell 
           * and 0 velocities. The number of electrons is equal to
           * the number of ions for hydrogen. 
           */ 
          x.set_e1( static_cast<Scalar>(j) + frand() );
          x.set_e2( static_cast<Scalar>(k) + frand() );
          
          // create electron macro particles
          pList.add(new Particle( x, u, eSpecies, TI_np2c, true ) );
          
          // create ion macro particles
          pList.add(new Particle( x, u, iSpecies, TI_np2c, true ) );
        }
        continue;
      }
      // If all neutral atoms have been ionized and we reach this point,
      // then we must create the final macro-particle for this cell.
      // Otherwise, we should continue with the for-loop.
      if ( (stmp*cellVol < 1.0) && (excessNumIons > 1.0) ) {
        // ionize the rest of the neutral atoms and place them in 
        // a single macro particle.
        x.set_e1( static_cast<Scalar>(j) + frand() );
        x.set_e2( static_cast<Scalar>(k) + frand() );
        
        // create electron macro particles
        pList.add(new Particle( x, u, eSpecies, excessNumIons, true ) );
        
        // create ion macro particles
        pList.add(new Particle( x, u, iSpecies, excessNumIons, true ) );

        stmp = 0.0;
        ptrNGD->set(j, k,  stmp);
        excessNumIons = 0;
        try{
          ptrNGD->setNGD_excessNumIons(j, k, excessNumIons);
        }
        catch(Oops& oops){
          oops.prepend("CsMCTI::tunnelIonizeCs: Error: \n"); //OK
          throw oops;
        }
      } 
    }
  }
  //
  // distribute the created particles during the tunneling ionization
  // in the appropriate particle group lists according to species
  // 
  region->addParticleList(pList);
}

/** 
 * A helper function used for the tunneling ionization of the 
 * Cs+ ions.
 */
void CsMCTI::tunnelIonizeCsPlus(ParticleGroupList& pgList, 
                                ParticleList& pList) 
  throw(Oops){
  //
  // check first if the ParticleGroupList is empty, if yes => return
  // 
  if (!pgList.head) 
    return;

  int n; // counter
  Scalar EfieldMagnitude;
  Scalar Probability;
  Scalar thisParticle_np2c;
  Vector2 thisParticle_x;
  Vector3 thisParticle_u;  
  /**
   * calculate the time step in atomic units of time
   */
  Scalar timeStep  = region->get_dt() / ATOMIC_TIME_UNIT;
  
  // 
  // loop over the particle group list
  // 
  oopicListIter<ParticleGroup> pgIter(pgList);
  for (pgIter.restart(); !pgIter.Done(); pgIter++) {
    //
    // loop over the particles in each particle group
    //
    int numParticles = pgIter()->get_n(); // number of particles in the group
    n = 0; 
    // 
    // need to keep track if the group becomes empty to
    // schedule its deletion
    // 
    BOOL IsParticleGroupEmpty; 
    while ( n < numParticles ) {
      thisParticle_x    = pgIter()->get_x(n);
      EfieldMagnitude = ((region->get_fields())->E(thisParticle_x)).magnitude()
        /ATOMIC_E_FIELD_UNIT_MKS;
      try{
        Probability = TIProbability(fp, EfieldMagnitude, timeStep, Zti[1]-1); 
      }
      catch(Oops& oops2){
        oops2.prepend("CsMCTI::tunnelIonizeCsPlus: Error: \n");//done
        throw oops2;
      }

      //
      // generate a "quasi" random number QRN and ionize the macroparticle if
      // QRN < Probability:
      // 
      if ( frand() < Probability ) {
        //
        // ionize the particle; note that in this case n is not incremented
        // since the last particle in the group is placed at the index n 
        // of the particle that we are going to ionize and delete from this
        // group
        // 
        thisParticle_np2c = pgIter()->get_np2c(n);
        thisParticle_u    = pgIter()->get_u(n);
        //
        // create electron macro particle
        //
        pList.add(new Particle( thisParticle_x, thisParticle_u, 
                                eSpecies, thisParticle_np2c, true ) );
        //
        // create CsPlusPlus ion macro particle
        // 
        pList.add(new Particle( thisParticle_x, thisParticle_u, 
                                iSpeciesPlusPlus, thisParticle_np2c, true ) );
        //
        // delete the old macro particle and keep track if the group is empty
        //
        IsParticleGroupEmpty = pgIter()->remove(n);
        numParticles--;
      } else {
        // 
        // this particle was not ionized => increment the "n" counter to
        // move to the next macro particle in the group
        // 
        n++;
      }
    }
  }
}

/**
 * a helper function to calculate the probability for tunneling 
 * ionization of Cs atom/ions from the ground state
 */
Scalar CsMCTI::TIProbability(Scalar (CsMCTI::*P)(Scalar&, Scalar&, int) const,
                            Scalar& E, Scalar& dt, int Zindex) 
  throw(Oops){
  Scalar Probability = 0.0;
  
  if ( E > EmaxTI[Zindex] ) { 
    //
    // if the field is greater than the max E field determined
    // from the inflection formula for the probability rate
    // then set the probability to 1. 
    // 
    Probability = 1.0;
  } else if ( E > EminTI[Zindex] ) {
    //
    // if E > EminTI (as determined from Keldish' condition)
    // and is less than EmaxTI then use the TI formula to
    // calculate the probability for ionization
    //
    try{
      Probability = (this->*P)(E, dt, Zindex);
    }
    catch(Oops& oops){
      oops.prepend("CsMCTI::TIProbability: Error:\n"); //CsMCTI::tunnelIonizeCsPlus
      throw oops;
    }

    if ( Probability > 1.0 ) 
      Probability = 1.0;
  } else {
    // 
    // this is the default case: for E < EminTI
    // the probability for TI is negligible, so set it to zero
    // 
    Probability = 0.0;
  }
  return Probability; 
}