limcti.cpp

pic 模拟程序！面向对象

      catch(Oops& oops2){
        oops2.prepend("LiMCTI::tunnelIonize: 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("LiMCTI::tunnelIonize: 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("LiMCTI::tunnelIonize: Error: \n");//MCTI::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("LiMCTI::tunnelIonize: 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("LiMCTI::tunnelIonize: 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 to calculate the probability for tunneling 
 * ionization of Hydrogen ions from the ground state
 */
Scalar LiMCTI::TIProbability(Scalar (LiMCTI::*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("LiMCTI::TIProbability: Error:\n");  //LiMCTI::tunnelIonize
      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; 
}