sptlrgn.cpp

pic 模拟程序！面向对象

   * 
   * D. A. Dimitrov 06/27/01. 
   */
  if (mcti) {
    try{
      mcti->tunnelIonize(particleGroupList);
    }
    catch(Oops& oops3){
      oops3.prepend("SpatialRegion::particleAdvance: Error: \n");
      throw oops3;
    }
  }

  // Need to pass particles across SpatialRegionBoundary's
  // This should go inside the loop when emit goes inside.
  //  ApplyToList(passParticles(),*boundaryList,Boundary);
		
  globalParticleLimit();
  nParticleFlag = FALSE;
  if (!fields->get_freezeFields()) {
    fields->compute_rho(particleGroupList,nSpecies);  //Get the charge density
    try{
      ApplyToList(putCharge_and_Current(t),*boundaryList,Boundary);
    }
    catch(Oops& oops2){
      oops2.prepend("SpatialRegion::particleAdvance: Error: \n");//OK
      throw oops2;
    }

    fields->ZeroCharge();  //not rho but charge from dielectrics, for diag only
    ApplyToList(putCharge(),*boundaryList,Boundary);   //for diag only
  }
}

void SpatialRegion::fieldsAdvance() {
  fields->advance(t, dt);//	advance fields
  t += dt*getFieldSubFlag();
  if (getBoltzmannFlag()) fields->IncBoltzParticles(
                                                    *particleGroupList[fields->get_bSpecies()->getID()]);
  fields->toNodes(t);
}


void SpatialRegion::emit(Species *species) 
  throw(Oops){
  oopicListIter<Boundary> nextb(*boundaryList);
  for (nextb.restart(); !nextb.Done(); nextb++){ // loop through boundaries
    try{
      addParticleList(nextb.current()->emit(t,getFieldSubFlag()*dt, species));
    }
    catch(Oops& oops){
      oops.prepend("SpatialRegion::emit: Error: \n"); //particleAdvance()
      throw oops;
    }
  }
}

void SpatialRegion::addParticleList(ParticleList& particleList)
{
  // empty Particle list into available ParticleGroups
  oopicListIter<Particle> pIter(particleList);
  for (pIter.restart(); !pIter.Done(); pIter++)
    {
      Particle* p = pIter();
      int s = p->get_speciesID();
      if (s < nSpecies)
        {
          oopicListIter<ParticleGroup> pgIter(*particleGroupList[s]);
          for (pgIter.restart(); !pgIter.Done(); pgIter++)
            if (pgIter()->add(p)) break;
          if (pgIter.Done())			//	no space; add new group
            {
              ParticleGroup* pg;
              if(getInfiniteBFlag() && getElectrostaticFlag())
                pg = new ParticleGroupIBES(maxN, p->get_speciesPtr(), p->get_np2c(), 
                                           p->isVariableWeight());
              else if (getInfiniteBFlag()) // relativistic Infinite B mover
                pg = new ParticleGroupIBEM(maxN, p->get_speciesPtr(), p->get_np2c(), 
                                           p->isVariableWeight());  
              else if (getElectrostaticFlag()) // add ES group
                pg = new ParticleGroupES(maxN, p->get_speciesPtr(), p->get_np2c(),
                                         p->isVariableWeight());
              else if (getSynchRadiationFlag()) // include synchrotron radiation
                pg = new ParticleGroupSR(maxN, p->get_speciesPtr(), p->get_np2c(), 
                                         p->isVariableWeight());
              else if(getNonRelFlag())
                pg = new ParticleGroupNR(maxN, p->get_speciesPtr(), p->get_np2c(), 
                                         p->isVariableWeight());
              else        //	add EM group
                pg = new ParticleGroup(maxN, p->get_speciesPtr(), p->get_np2c(), 
                                       p->isVariableWeight());
              pg->add(p);
              particleGroupList[s]->add(pg);
            }
        }
      pIter.deleteCurrent(); // remove list element and data
    }





}

//---------------------------------------------------------------------
// packParticleGroups() repacks the groups to minimize the number of groups.
// It will result in n full groups plus a partially full remainder group,
// which gets placed at the head of the list. CURRENTLY NON-FUNCTIONAL

void SpatialRegion::packParticleGroups(int i)
{
//  ParticleGroup *current;
  oopicListIter<ParticleGroup> source(*particleGroupList[i]);
  oopicListIter<ParticleGroup> dest(*particleGroupList[i]);
  for (source.restart(); !source.Done(); source++){
    if (!source()->vary_np2c() && source()->get_n() < source()->get_maxN())
      for (dest.restart(); !dest.Done(); dest++){
        if (source() != dest()) 
          if (dest()->add(*source())) break; // TRUE = finished adding
      }
    if (!source()->get_n()) source.deleteCurrent();
  }
}

void SpatialRegion::increaseParticleWeight(void)
{
  for (int i=0; i<nSpecies; i++)
    {
      oopicListIter<ParticleGroup> nextpg(*particleGroupList[i]);
      for (nextpg.restart(); !nextpg.Done(); nextpg++)
        nextpg()->increaseParticleWeight();
    }
   oopicListIter<Boundary>	nextb(*boundaryList);
  for(nextb.restart(); !nextb.Done(); nextb++)
    nextb()->increaseParticleWeight();
}

void SpatialRegion::increaseParticleWeight(int i)
{

  oopicListIter<ParticleGroup> nextpg(*particleGroupList[i]);
  for (nextpg.restart(); !nextpg.Done(); nextpg++)
    nextpg()->increaseParticleWeight();

  //	packParticleGroups(i);

  oopicListIter<Boundary>	nextb(*boundaryList);
  for(nextb.restart(); !nextb.Done(); nextb++)
    nextb()->increaseParticleWeight(i);
}



#ifdef HAVE_HDF5
int SpatialRegion::dumpH5(dumpHDF5 &dumpObj)
{
  std::cout << "Entering spatialRegion::DumpH5"<<endl;

  static int end=-1;
//  int nn;
//  FILE *DropFile;
//  char dropfile[80];
  int seqNumber = 0;
  int mpiRank = 1;
  hid_t fileId, groupId;// datasetId,dataspaceId;
  herr_t status;
  hsize_t rank; 
  hsize_t     dimsf[1];   
  hsize_t dims;
  hid_t scalarType = dumpObj.getHDF5ScalarType();
  hid_t  attrdataspaceId,attributeId; //attrDS;
  int nnnn;

 

  fileId = H5Fopen(dumpObj.getDumpFileName().c_str(), H5F_ACC_RDWR, H5P_DEFAULT);  
  hid_t bGroupId = H5Gcreate(fileId,"/Boundaries",0);




  groupId = H5Gcreate(fileId,"/SpatialRegion",0);
//     oopicListIter<Diag>	nextd(*diagList);
//    nnnn = 0;
//    for (nextd.restart(); !nextd.Done(); nextd++)
//      {
//        if (!(nextd.current()->dumpH5(dumpObj,nnnn)))
//  	cout << "Dump: dump failed to save diagnostics." << endl;
//        nnnn++;
//      }
  rank = 1;
  dimsf[0] = 4;
  string datasetName = "dims";  



// dump dimensions of global grid
 //   dims = 4;
//  #ifdef MPI_VERSION
//    int globalGridSize[4] = {grid->MaxMPIJ,grid->MaxMPIK,grid->M,grid->N};
//  #else
//    int globalGridSize[4] = {grid->getJ(),grid->getK(),1,1};
//  #endif 
//    attrdataspaceId = H5Screate_simple(1, &dims, NULL);
//    attributeId = H5Acreate(groupId, "globalGridSize",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
//    status = H5Awrite(attributeId, H5T_NATIVE_INT, globalGridSize);
//    status = H5Aclose(attributeId);
//    status = H5Sclose(attrdataspaceId);
//    cerr << "dumped global grid dimensions: "<< globalGridSize[0]<<"  "<< globalGridSize[1]<<endl;
   
//  // dump positions of grid
//    int J = grid->getJ();
//    int K = grid->getK();
//    Vector2 **X = grid->getX();
//    Scalar *x1List = new Scalar[J+1];
//    Scalar *x2List = new Scalar[K+1];
  
//    for(int i=0;i<=J;i++)
//      x1List[i] = X[i][0].e1();
  
//    for(int i=0;i<=K;i++)
//      x2List[i] = X[0][i].e2();

//    dimsf[0] = J+1;
//    dataspaceId = H5Screate_simple(rank, dimsf, NULL); 
//    datasetId = H5Dcreate(groupId, "x1Nodes", scalarType, dataspaceId, H5P_DEFAULT);
//    status = H5Dwrite(datasetId, scalarType, H5S_ALL,H5S_ALL , H5P_DEFAULT, x1List);
//    status = H5Sclose(dataspaceId); 
//    status = H5Dclose(datasetId);
  
//    dimsf[0] = K+1;
//    dataspaceId = H5Screate_simple(rank, dimsf, NULL); 
//    datasetId = H5Dcreate(groupId, "x2Nodes", scalarType, dataspaceId, H5P_DEFAULT);
//    status = H5Dwrite(datasetId, scalarType, H5S_ALL,H5S_ALL , H5P_DEFAULT, x2List);
//    status = H5Sclose(dataspaceId); 
//    status = H5Dclose(datasetId);

//    delete[] x1List;
//    delete[] x2List;

// dump random seed  as attribute
  dims = 1;
  attrdataspaceId = H5Screate(H5S_SCALAR);
  attributeId = H5Acreate(groupId, "frandseed",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT, &frandseed);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);
  // cerr << "dumped random seed: "<<frandseed<<endl;
// dump time as attribute
  dims = 1;
  attrdataspaceId = H5Screate(H5S_SCALAR);
  attributeId = H5Acreate(groupId, "simulationTime",scalarType, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, scalarType, &t);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);
  
  int temp = boundaryList->nItems();

// dump number of boundaries as attribute
  dims = 1;
//attrdataspaceId = H5Screate_simple(1, &dims, NULL);
  
  attrdataspaceId = H5Screate(H5S_SCALAR);
  attributeId = H5Acreate(groupId, "numBoundaries",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT, &temp );
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);

  temp = diagList->nItems();

// dump number of diagnostics as attribute
  dims = 1;
  attrdataspaceId = H5Screate(H5S_SCALAR);
  attributeId = H5Acreate(groupId, "numDiagnostics",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT, &temp  );
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);

// dump number of species  as attribute
  dims = 1;
  attrdataspaceId = H5Screate(H5S_SCALAR);
  attributeId = H5Acreate(groupId, "nSpecies",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
  status = H5Awrite(attributeId, H5T_NATIVE_INT, &nSpecies);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);

// dump dimensions of region as attribute
  dims = 4;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);

  attributeId = H5Acreate(groupId, "dims",scalarType, attrdataspaceId, H5P_DEFAULT);
   
  Scalar xw[4];  // in this format:  xs, ys, xf, yf
  Vector2 **X=grid->getX();
  xw[0] = X[0][0].e1();
  xw[1] = X[0][0].e2();
  xw[2] = X[getJ()][getK()].e1();
  xw[3] = X[getJ()][getK()].e2();
  // status = H5Dwrite(datasetId, scalarType, H5S_ALL, H5S_ALL,
  //	      H5P_DEFAULT, xw);
  
  status = H5Awrite(attributeId, scalarType, &xw);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);



  H5Gset_comment(fileId,"/Boundaries","test comment");







//  // dump dimensions of region as attribute
//    dims = 3;
//    attrdataspaceId = H5Screate_simple(1, &dims, NULL);

//    attributeId = H5Acreate(groupId, "dims",H5T_NATIVE_STRING, attrdataspaceId, H5P_DEFAULT);
   
//    Scalar xw[4];  // in this format:  xs, ys, xf, yf
//    Vector2 **X=grid->getX();
//    xw[0] = X[0][0].e1();
//    xw[1] = X[0][0].e2();
//    xw[2] = X[getJ()][getK()].e1();
//    xw[3] = X[getJ()][getK()].e2();
//    // status = H5Dwrite(datasetId, scalarType, H5S_ALL, H5S_ALL,
//    //	      H5P_DEFAULT, xw);
  
//    status = H5Awrite(attributeId, scalarType, &xw);
//    status = H5Aclose(attributeId);
//    status = H5Sclose(attrdataspaceId);


  hid_t       file_id; //group1_id, group2_id, group3_id;  /* identifiers */
  hid_t       aid, atype, attr;
   
  
  /* Create a new file using default properties. */
  file_id = fileId;
  
  dims = 4;
  aid = H5Screate(H5S_SCALAR);
  atype = H5Tcopy (H5T_C_S1);
  
  status = H5Tset_size (atype, 80*nSpecies);
  
  attr = H5Acreate (bGroupId, "data_contents", atype, aid, H5P_DEFAULT);
  
  char *what = new char[80*nSpecies];
  if (nSpecies!=0){
    strcpy(what, "");
  }
    
  // write a comma delimited list of the species names into an attribute
   oopicListIter<Species> spIter(*speciesList);
   spIter.restart();
   if(!spIter.Done())
     {
       strcat(what,spIter.current()->get_name().c_str() );
       spIter++;
     }

   while(!spIter.Done())
     {
       strcat(what,",");
       strcat(what,spIter.current()->get_name().c_str() );
       spIter++;
       cout<<"count"<<what<<endl;
     }
   cout<<what<<endl;
  

   
//     strstream spec;
//      spec << "/" << speciesIter.current()->get_name().c_str() <<ends;    
//      string spGroupName = spec.str();



   status = H5Awrite (attr, atype,what );

   //status = H5Gclose(bGroupId);

   status = H5Sclose (aid);

   status = H5Aclose (attr);

   delete [] what;
  //    // make a drop file entry
  //    nn = strlen(theDumpFile);               
  //    strcpy(dropfile, theDumpFile);
  //    dropfile[nn-2] = 'r';
  //    if ((DropFile = fopen(dropfile, "a")) == NULL)
  //      {
  //        printf("filename = %s\n", dropfile);
  //        puts("Drop: open failed error#1");
  //        return (0);
  //      }
  //    fprintf(DropFile, "%g ", t);
  
  //    for (int i=0; i<nSpecies; i++){
  //      int number = CountParticles(i);
  //      fprintf(DropFile, " %d ", number);
  //    }
  
  //    // Drop the number of shifts the spatial region has done.
  //    fprintf(DropFile," %lu ", shiftNum);
  
  //    fprintf(DropFile, "\n");
  //    fclose(DropFile);
 
//  Dump Boundaries
  oopicListIter<Boundary>	nextb(*boundaryList);
  nnnn = 0;
  for (nextb.restart(); !nextb.Done(); nextb++)
    {
      // string groupName = nn;
      // hid_t groupId = H5Gcreate(fileId, groupName.c_str(),0);
      if (!(nextb.current()->dumpH5(dumpObj,nnnn)))
        printf("Dump: dump failed to save boundaries");
      nnnn++;
    }
	
 H5Gclose(bGroupId);

//  Dump diagnostics
  //   {
  //      oopicListIter <Diagnostics> Diter(*theDiagnostics);
  //      for(Diter.restart();!Diter.Done();Diter++)
  //        Diter.current()->dumpH5(dumpObj);
  //    }  

//
// dump fields
// 

 //  int J = grid->getJ();
//    int K = grid->getK();

//   for(int j=0;j<J;j++){
//      for(int k=0;k<K;k++){
//        cout << fields->getIntEdl()[j][k].e1()<<" ";
//        //  cout << fields->getENode()[j][k].e1()<<"\n";
//      }
//    }
 
 fields->dumpH5(dumpObj);
//
// dump NGD
// 
  oopicListIter<NGD> NGDIter(NGDList);
  nnnn = 0;
  for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ){ 
    NGDIter.current()->dumpH5(dumpObj,nnnn);
    nnnn++;
  }

// 
// dump Species
// 

  int isp= 0;
  oopicListIter<Species> speciesIter(*speciesList);
  for (speciesIter.restart(); !speciesIter.Done(); speciesIter++){
        
    int specID = speciesIter.current()->getID();
    int ptclGrpCount = 0;

    //    strstream sp;
#ifdef HAVE_SSTREAM
    //stringstream sp;
#else
    //strstream sp;
#endif

    //sp << "/" << speciesIter.current()->get_name().c_str() <<ends;    
    //string spGroupName = sp.str();
    //   hid_t spGroupId = H5Gcreate(fileId,spGroupName.c_str(),0);
     hid_t spGroupId = H5Gcreate(fileId,speciesIter.current()->get_name().c_str(),0);

     oopicListIter<ParticleGroup> pgscan(*particleGroupList[specID]);
     for(pgscan.restart();!pgscan.Done();pgscan++){