sptlrgn.cpp

pic 模拟程序！面向对象

       status = H5Aclose(attrId);

      for(int j=0;j<nPtclGrp;j++) 
	{
          //	  strstream s;
#ifdef HAVE_SSTREAM
         stringstream s;
#else
         strstream s;
#endif

	  s << spGrpName <<"/pGrp"<<j<<ends;  
	  string pGrpName = s.str();
	  hid_t pGroupId = H5Gopen(fileId,pGrpName.c_str());
// Read the speciesID attribute
	  attrId = H5Aopen_name(pGroupId, "speciesID");
	  status = H5Aread(attrId,H5T_NATIVE_INT , &speciesID);
	  status = H5Aclose(attrId);
      
//  	  oopicListIter<Species> speciesIter(*speciesList);
//  	  for (speciesIter.restart(); !speciesIter.Done(); speciesIter++)
//  	    if (speciesIter.current()->getID() == speciesID)
//  	      species = speciesIter.current();
      
// Read the vary_np2c attribute
	  attrId = H5Aopen_name(pGroupId, "vary_np2c");
	  status = H5Aread(attrId,H5T_NATIVE_INT , &vary_np2c);
	  status = H5Aclose(attrId);
      
// Read the np2c attribute
	  attrId = H5Aopen_name(pGroupId, "np2c");
	  status = H5Aread(attrId,scalarType , &np2c0);
	  status = H5Aclose(attrId);

// Making a new particleGroup
	  np2c0 *= np2cFactor/(Scalar(1+duplicateParticles));;
	  ParticleGroup* pg;
	  if(getInfiniteBFlag() && getElectrostaticFlag())
	    pg = new ParticleGroupIBES(maxN, species, np2c0, vary_np2c);
	  else if (getInfiniteBFlag()) // relativistic Infinite B mover
	    pg = new ParticleGroupIBEM(maxN, species, np2c0, vary_np2c);  
	  else if (getElectrostaticFlag()) // add ES group
	    pg = new ParticleGroupES(maxN, species, np2c0, vary_np2c);
	  else if (getSynchRadiationFlag()) // add SR group
	    pg = new ParticleGroupSR(maxN, species, np2c0, vary_np2c);
	  else if(getNonRelFlag())
	    pg = new ParticleGroupNR(maxN, species, np2c0, vary_np2c);
	  else        //	add EM group
	    pg = new ParticleGroup(maxN, species, np2c0, vary_np2c);

	  pg->restoreH5(restoreObj, grid, np2cFactor,pGrpName);
	  
	  ParticleGroup** k;
	  if(duplicateParticles > 0) 
	    { 
	      fprintf(stderr, "\nduplicateParticles = %d", duplicateParticles);
	      k = pg->duplicateParticles(duplicateParticles); 
	      if(k != NULL)
		{
		  for(int i=0; i < duplicateParticles; i++)
		    {
		      particleGroupList[speciesID]->add(k[i]);
		    }
		}
	      else
		{
		  fprintf(stderr, "\nError duplicating particles");
		}
	    }
	  
	  particleGroupList[speciesID]->add(pg);
	  H5Gclose(pGroupId);
	}
      H5Gclose(spGroupId);
    }
 
  H5Gclose(groupId);
  H5Fclose(fileId);
  return(1);
}
#endif

int SpatialRegion::Restore(FILE* DMPFile)
{
  Scalar x1min, x2min, x1max, x2max;
  Scalar gx1min, gx2min, gx1max, gx2max;
  Vector2 **X=grid->getX();

  int N_Boundaries;
//  int found;
  int i;

  XGRead(&frandseed, 4, 1, DMPFile, "long");
	
  int ReadNameLength;
  char buf[512];
	
  XGRead(&ReadNameLength,4,1,DMPFile,"int");   // Read length of the region name
  XGRead(buf,1,ReadNameLength+1,DMPFile,"char"); // Read region name

  XGRead(&t, ScalarInt, 1, DMPFile, ScalarChar); // Read simulation time

  //  Restore Boundaries must be in the same order as dumped
  // Should give them titles as in diagnostics


  // read in the dimensions of the dumpfile being read
  XGRead(&x1min, ScalarInt, 1, DMPFile, ScalarChar);
  XGRead(&x2min, ScalarInt, 1, DMPFile, ScalarChar);
  XGRead(&x1max, ScalarInt, 1, DMPFile, ScalarChar);
  XGRead(&x2max, ScalarInt, 1, DMPFile, ScalarChar);

  // decide if this is in our domain so we need to restore some of it
  gx1min = X[0][0].e1();
  gx2min = X[0][0].e2();
  gx1max = X[getJ()][getK()].e1();
  gx2max = X[getJ()][getK()].e2();

  // all these checks exit if the file we're reading is out
  // of our domain.
  if(x1min > gx1max) return 1;
  if(x2min > gx2max) return 1;
  if(x1max < gx1min) return 1;
  if(x2max < gx2min) return 1;

	
  //read in the number of boundaries
  XGRead(&N_Boundaries, 4, 1, DMPFile, "int"); 

  // read in the dump information for each boundary
  for(i=0;i<N_Boundaries;i++) {
    Scalar A1,A2,B1,B2;  // phys. dimensions of boundary
    int BoundaryType; // Type of this boundary
    int WasRestored=0; //Did one of the boundaries restore from this boundary?
    int nQuads;  // number of ordered quads of floats to restore
		
    // read in the boundary extent
    XGRead(&A1,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&A2,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&B1,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&B2,ScalarInt,1,DMPFile,ScalarChar);

    XGRead(&BoundaryType,4,1,DMPFile,"int");
    // This next int tells us if there is anything to restore
    XGRead(&nQuads,4,1,DMPFile,"int");
    if(nQuads>0) {
      // try to find a matching boundary
      oopicListIter<Boundary>	nextb(*boundaryList);
      for (nextb.restart(); !nextb.Done()&& WasRestored==0; nextb++){
        WasRestored=nextb.current()->Restore(DMPFile,BoundaryType,A1,A2,B1,B2,nQuads);
				
      }
    }
    if(WasRestored==0)  { //nothing to restore:  load in the dummy stuff
      XGRead(&nQuads,4,1,DMPFile,"int");
      XGRead(&nQuads,4,1,DMPFile,"int");
      if(nQuads > 0) {  //wasn't restored, but there's data we have to read
        int l;
        Scalar QuadIn[4];
        for(l=0;l<nQuads;l++)
          XGRead(&QuadIn,ScalarInt,4,DMPFile,ScalarChar);
      }
    }
  }
	
  //  Restore diagnostic

  int N_Diags=0;
  XGRead(&N_Diags, 4, 1, DMPFile, "int");

  for(i=0;i<N_Diags;i++) {
    int WasRestored=0;
    Scalar A1,A2,B1,B2;
    int xlength;
    int ylength;
    int tlength;
    char title[256];
    int HistFlag;

    XGRead(&A1,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&A2,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&B1,ScalarInt,1,DMPFile,ScalarChar);
    XGRead(&B2,ScalarInt,1,DMPFile,ScalarChar);
		
    XGRead(&xlength,4,1,DMPFile,"int");
    XGRead(&ylength,4,1,DMPFile,"int");

    XGRead(&tlength,4,1,DMPFile,"int");

    XGRead(title,1,tlength+1,DMPFile,"char"); // +1 to get the null
    XGRead(&HistFlag,4,1,DMPFile,"int");

    if(HistFlag) {
      oopicListIter<Diag>	nextd(*diagList);
      for (nextd.restart(); !nextd.Done() && WasRestored==0; nextd++){
        WasRestored=nextd.current()->Restore(DMPFile,A1,A2,B1,B2,xlength,ylength,title);
      }
      if(WasRestored==0) { //nothing to restore:  load in the data
        int datlen;
        int dummy[11];
        int l;
				// first read the generic history stuff
        XGRead(&datlen,4,1,DMPFile,"int");
				// Read in two unused ints (hist_hi and alloc_size)
        XGRead(dummy,4,2,DMPFile,"int");
				// read in the unused data
        for(l=0;l<datlen;l++) { 
          Scalar dummy2;
          XGRead(&dummy2,ScalarInt,1,DMPFile,ScalarChar);
        }
				// Now read the particular history stuff
        XGRead(&datlen,4,1,DMPFile,"int");
				// Read in 11 unused ints 
        XGRead(dummy,4,11,DMPFile,"int");
				// read in the unused data
        for(l=0;l<datlen;l++) { 
          Scalar dummy2;
          XGRead(&dummy2,ScalarInt,1,DMPFile,ScalarChar);
        }
				


      }


    }
  }
 
  fields->Restore(DMPFile);
  //
  // restore of the neutral gas densities
  // 
  oopicListIter<NGD> NGDIter(NGDList);
  for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ) 
    NGDIter.current()->Restore(DMPFile); 

  int speciesID;
  Scalar np2c0;
  int vary_np2c;
  Species* species=0;

  XGRead(&nSpecies, 4, 1, DMPFile, "int");
  XGRead(&speciesID, 4, 1, DMPFile, "int");
  while (!(speciesID==-1)){

    oopicListIter<Species> speciesIter(*speciesList);
    for (speciesIter.restart(); !speciesIter.Done(); speciesIter++)
      if (speciesIter.current()->getID() == speciesID)
        species = speciesIter.current();
    
    XGRead(&vary_np2c, 4, 1, DMPFile, "int");
    XGRead(&np2c0, ScalarInt, 1, DMPFile, ScalarChar);
		
    // Making a new particleGroup
		
    np2c0 *= np2cFactor/(Scalar(1+duplicateParticles));
    ParticleGroup* pg;
    if(getInfiniteBFlag() && getElectrostaticFlag())
      pg = new ParticleGroupIBES(maxN, species, np2c0, vary_np2c);
    else if (getInfiniteBFlag()) // relativistic Infinite B mover
      pg = new ParticleGroupIBEM(maxN, species, np2c0, vary_np2c);  
    else if (getElectrostaticFlag()) // add ES group
      pg = new ParticleGroupES(maxN, species, np2c0, vary_np2c);
    else if (getSynchRadiationFlag()) // add SR group
      pg = new ParticleGroupSR(maxN, species, np2c0, vary_np2c);
    else if(getNonRelFlag())
      pg = new ParticleGroupNR(maxN, species, np2c0, vary_np2c);
    else        //	add EM group
      pg = new ParticleGroup(maxN, species, np2c0, vary_np2c);

    pg->Restore(DMPFile, grid, np2cFactor);
    ParticleGroup** k;
    if(duplicateParticles > 0) 
      { 
	//	fprintf(stderr, "\nduplicateParticles = %d", duplicateParticles);
	k = pg->duplicateParticles(duplicateParticles); 
	if(k != NULL)
	  {
	    for(int i=0; i < duplicateParticles; i++)
	      {
		particleGroupList[speciesID]->add(k[i]);
	      }
	  }
	else
	  {
	    fprintf(stderr, "\nError duplicating particles");
	  }
      }
    particleGroupList[speciesID]->add(pg);

    XGRead(&speciesID, 4, 1, DMPFile, "int");
  }

  return (1);
}

void SpatialRegion::shiftRegion() 
  throw(Oops){

#ifdef MPI_VERSION
  extern int MPI_RANK;
#endif

  //
  // Update the next shift time
  //
  next_shift_time += shift_dt;

  // Increment the shift count
  shiftNum++;

  //
  // Shift all of the species and send out
  //
  // cout << "MPI_RANK = " << MPI_RANK << ": shifting particles." << endl;
  oopicListIter<Species> siter(*speciesList);
  for (siter.restart(); !siter.Done(); siter++){

    // Get the species
    Species* species = siter.current();
    int i = species->getID();

    // Shift all of the particle groups for that species
    // This collects all of the particles into the boundaries
    oopicListIter<ParticleGroup> nextpg(*particleGroupList[i]);
    for (nextpg.restart(); !nextpg.Done(); nextpg++) 
      nextpg()->shiftParticles(*fields, shift);
  }

  //
  // Pass the particles to the other region.
  //

  // Determine the boundary that is shifted out
  Boundary* bPtrOut = 0;
  switch(shiftDir){
   
  case 1:
    bPtrOut=grid->GetBC2()[0][1];
    break;

  case 2:
    bPtrOut=grid->GetBC2()[grid->getJ()][1];
    break;

  case 3:
    bPtrOut=grid->GetBC1()[1][0];
    break;

  case 4:
    bPtrOut=grid->GetBC1()[1][grid->getK()];
    break;

  }

  // Pass the particles.  This does nothing except for spatial
  // region boundaries.  Presumably, none of the other regions
  // collected particles.
  if(!bPtrOut){
    static bool outPrinted = false;
    if(!outPrinted) cout << "Error - no outward boundary for shiftDir = " << 
                      shiftDir << endl;
    outPrinted = true;
  }
  else{
    // cout << "MPI_RANK = " << MPI_RANK << ": passing shifted particles." << endl;
    bPtrOut->passShiftedParticles();
  }

  // Determine the boundary that is shifted in.
  // Same as above, but switch ends.
  Boundary* bPtrIn = 0;
  switch(shiftDir){

  case 2:
    bPtrIn=grid->GetBC2()[0][1];
    break;

  case 1:
    bPtrIn=grid->GetBC2()[grid->getJ()][1];
    break;

  case 4:
    bPtrIn=grid->GetBC1()[1][0];
    break;

  case 3:
    bPtrIn=grid->GetBC1()[1][grid->getK()];
    break;

  }

  // Better have a boundary
  if(!bPtrIn){
    static bool inPrinted = false;
    if(!inPrinted) cout << "Error - no inward boundary for shiftDir = " << 
                     shiftDir << endl;
    inPrinted = true;
  }
  else{
    // cout << "MPI_RANK = " << MPI_RANK << ": asking for stripe." << endl;
    bPtrIn->askStripe();

    //
    // Send fields to right, then shift internal
    // cout << "MPI_RANK = " << MPI_RANK << ": sending fields." << endl;
    //
    // We want to shift the data associated with the neutral gas density
    // when we do the Fields send and shift so we will pass a pointer to 
    // the neutral gas density. dad 05/16/01.
    //
 
    if(bPtrOut) fields->sendFieldsAndShift(shiftDir, bPtrOut);

    //
    // Receive any particles from the boundary
    //
    if(bPtrIn->getBCType() == SPATIAL_REGION_BOUNDARY){
      // cout << "MPI_RANK = " << MPI_RANK << ": receiving particles." << endl;
      ParticleList& pl = bPtrIn->recvShiftedParticles();
      // cout << "MPI_RANK = " << MPI_RANK << ": adding particles." << endl;
      addParticleList(pl);	// This empties the particle list
    }
    
    /**
     * We want to follow the approach taken in communicating the 
     * fields when transfering the NGD's. However, we will do this
     * work only if the oopicList with NGD's is not empty. dad-05.21.01.
     */
    if ( !NGDList.isEmpty() ) {
      /**
       * ask first for a stripe of NGDs. 
       */

      bPtrIn->askNGDStripe(); 

      //
      // Send NGDs to right, then shift internal. 
      //
      if(bPtrOut) 
        try{
          ptrMapNGDs->sendNGDsAndShift(shiftDir, bPtrOut);
        }
        catch(Oops& oops2){
          oops2.prepend("SpatialRegion::shiftRegion: Error: \n");//done
          throw oops2;
        }
     }
  }

  //
  // Load any shiftLoads
  //
  // cout << "MPI_RANK = " << MPI_RANK << ": loading shift loads." << endl;
  oopicListIter<Load> nextload(*theLoads);
  for (nextload.restart(); !nextload.Done(); nextload++)
    if(nextload()->name == ostring("shiftLoad")) {
      try{
        nextload()->load_it();
      }
      catch(Oops& oops){
        oops.prepend("SpatialRegion::shiftRegion: Error: \n"); //SpatialRegion::advance
        throw oops;
      }
    }
  //
  // Collect the shifted fields
  //
  // cout << "MPI_RANK = " << MPI_RANK << ": receiving stripe." << endl;
  if(bPtrIn) fields->recvShiftedFields(shiftDir, bPtrIn,t,dt);
  // cout << "MPI_RANK = " << MPI_RANK << ": shift completed." << endl;

  /**
   * Collect the shifted NGDs
   */
  if ( !NGDList.isEmpty() ) {
    if ( bPtrIn )
      ptrMapNGDs->recvShiftedNGDs(shiftDir, bPtrIn);
  }

}