fields.cpp

pic 模拟程序！面向对象

      theBoltzmann->Restore(DMPFile);

   return(1);
}


#ifdef HAVE_HDF5
int Fields::restoreH5(dumpHDF5 &restoreObj)
{

  //cerr << "entering Fields::restoreH5()\n";

#ifdef MPI_VERSION
  extern int MPI_RANK;
  //  cerr << " ,MPI_RANK "<<MPI_RANK<<"   "<<restoreObj.getDumpFileName()<<"\n";
  
#endif
  
  hid_t fileId, groupId, datasetId,dataspaceId;
  herr_t status;
  string outFile;
//  hsize_t rank; 
//  hsize_t     dimsf[1];   
//  hsize_t dims;
  hid_t scalarType = restoreObj.getHDF5ScalarType();
//  hid_t  attrspaceId,
  hid_t attrId;
//  hid_t bMetaGroupId; 
//  hssize_t offset[2]; /* Offset of hyperslab */
//  hsize_t count[2];  /* Block count */
//  hssize_t offsetFile[3]; /* Offset of hyperslab */
//  hsize_t countFile[3];  /* Block count */
//  int EMdampingFlag;
  int jl,kl,Jl,Kl;
  //Scalar data,dataE,dataB,dataI;

  fileId = H5Fopen(restoreObj.getDumpFileName().c_str(), H5F_ACC_RDONLY, H5P_DEFAULT);

  groupId = H5Gopen(fileId, "/Fields");
  
  
  if(!ElectrostaticFlag) {
      
//get size of Fields:
    datasetId = H5Dopen(groupId,"positions");
    dataspaceId = H5Dget_space(datasetId );
    int dataRank = H5Sget_simple_extent_ndims(dataspaceId);
    hsize_t *sizes = new hsize_t[dataRank];
    int res = H5Sget_simple_extent_dims(dataspaceId, sizes, NULL);
    Jl = sizes[0] -1;
    Kl = sizes[1] - 1; //the -1 is to set the max index
    
    int third  = sizes[2];
    
    
    hsize_t dim2 = 2;
    // get dataspaces fro fields
    
    
    hid_t EdlsetId = H5Dopen(groupId, "Edl");
    hid_t EdlspaceId = H5Dget_space(EdlsetId);
    dataRank = H5Sget_simple_extent_ndims(EdlspaceId);
    delete [] sizes;
    sizes = new hsize_t[dataRank];
    res = H5Sget_simple_extent_dims(EdlspaceId, sizes, NULL);
    
    hid_t BdSsetId = H5Dopen(groupId, "BdS");
    hid_t BdSspaceId = H5Dget_space(BdSsetId);
    
    hid_t IsetId = H5Dopen(groupId, "I");
    hid_t IspaceId = H5Dget_space(IsetId);
    
    
    Scalar *pos = new Scalar[(Jl+1)*(Kl+1)*2];
    status = H5Dread(datasetId, scalarType,H5S_ALL, dataspaceId, H5P_DEFAULT, pos);
    
    
    Scalar *data = new Scalar[(Jl+1)*(Kl+1)*3];
    status = H5Dread(EdlsetId, scalarType,H5S_ALL, EdlspaceId, H5P_DEFAULT, data);
    
    Scalar *BdSdata = new Scalar[(Jl+1)*(Kl+1)*3];
    status = H5Dread(BdSsetId, scalarType,H5S_ALL, BdSspaceId, H5P_DEFAULT, BdSdata);
    
    Scalar *Idata = new Scalar[(Jl+1)*(Kl+1)*3];
    status = H5Dread(IsetId, scalarType,H5S_ALL, IspaceId, H5P_DEFAULT, Idata);
    
    
    //ofstream out;
  int compNum = 3;
  int maxNum = (Jl+1)*(Kl+1);
  
  for (jl=0; jl<=Jl;jl++){
    for (kl=0; kl<=Kl;kl++)
      {
	Scalar A1,A2;  // coordinates of the fields read in
	Vector2 Xloc;
	Vector3 E,B,EN,BN,Il;
	
	A1 = pos[(Kl+1)*2*jl + kl*2];
	A2 = pos[(Kl+1)*2*jl + kl*2 + 1];
	Xloc = grid->getNearestGridCoords(Vector2(A1,A2));
	
	
	
	if(Xloc.e1() >= 0) // valid location
	  {
	    int ji,ki;
	    // the following operations make sure
	    // stuff doesn't get rounded down when
	    // converted from float to int
	    ji = (int)(Xloc.e1()+0.5);
	    ki = (int)(Xloc.e2()+0.5);
	    
	    intEdl[ji][ki].set_e1(data[(Kl+1)*3*jl + kl*3]);
	    intEdl[ji][ki].set_e2(data[(Kl+1)*3*jl + kl*3 + 1]);
	    intEdl[ji][ki].set_e3(data[(Kl+1)*3*jl + kl*3 + 2]);
	    
	    intBdS[ji][ki].set_e1(BdSdata[(Kl+1)*3*jl + kl*3]);
	    intBdS[ji][ki].set_e2(BdSdata[(Kl+1)*3*jl + kl*3 + 1]);
	    intBdS[ji][ki].set_e3(BdSdata[(Kl+1)*3*jl + kl*3 + 2]);
	    
	    I[ji][ki].set_e1(Idata[(Kl+1)*3*jl + kl*3]);
	    I[ji][ki].set_e2(Idata[(Kl+1)*3*jl + kl*3 + 1]);
	    I[ji][ki].set_e3(Idata[(Kl+1)*3*jl + kl*3 + 2]);
	  }
      }
  }
  
  
      delete[] BdSdata;
      delete[] Idata;

      status = H5Sclose(BdSspaceId);
      status = H5Dclose(BdSsetId);
      status = H5Sclose(IspaceId);
      status = H5Dclose(IsetId);
      status =    H5Sclose(dataspaceId);
      status =    H5Dclose(datasetId);
      status = H5Sclose(EdlspaceId);
      status = H5Dclose(EdlsetId);

      // Read the EM damping flag attribute

      attrId = H5Aopen_name(groupId, "EMdampingFlag");
      Scalar EMdampingFlag;
      status = H5Aread(attrId,scalarType ,&EMdampingFlag);
      status = H5Aclose(attrId);
    
      if(EMdampingFlag>0){

	EdlsetId = H5Dopen(groupId,"EdlBar" );
	EdlspaceId = H5Dget_space(EdlsetId );
	status = H5Dread(EdlsetId, scalarType,H5S_ALL, EdlspaceId, H5P_DEFAULT, data);


// read in the EMDamping stuff
        for (jl=0; jl<=Jl;jl++)
          for (kl=0; kl<=Kl;kl++)
            {
              Scalar A1,A2;  // coordinates of the fields read in
              Vector2 Xloc;
              Vector3 Emd;
	      A1 = pos[(Kl+1)*2*jl + kl*2];
	      A2 = pos[(Kl+1)*2*jl + kl*2 + 1];

	      // find the grid locations here
              Xloc = grid->getNearestGridCoords(Vector2(A1,A2));
              //Emd
              if(Xloc.e1() >= 0) // valid location
                {
                  int ji,ki;
                  // the following operations make sure
                  // stuff doesn't get rounded down when
                  // converted from float to int
                  ji = (int)(Xloc.e1()+0.5);
                  ki = (int)(Xloc.e2()+0.5);
		 
		  intEdlBar[ji][ki].set_e1(data[(Kl+1)*3*jl + kl*3]);
		  intEdlBar[ji][ki].set_e2(data[(Kl+1)*3*jl + kl*3 + 1]);
		  intEdlBar[ji][ki].set_e3(data[(Kl+1)*3*jl + kl*3 + 2]);
		}
            }			
	
// if we've just turned EMDamping on.
        if (EMdamping != EMdampingFlag){
          cout << "**********************warning**********************" << endl;
          cout << "EMdamping value has been turned on, using damping." << endl;
          cout << "Weighted average is beening started with the instantaneous" << endl;
          cout << "value of the fields." << endl;
          for (int j=0; j<=J;j++)
            for (int k=0; k<=K;k++)
              intEdlBar[j][k] = intEdl[j][k];
        }
	
	status = H5Sclose(EdlspaceId);
	status = H5Dclose(EdlsetId);
	
      }
      delete [] sizes;
  }
  H5Gclose(groupId);
  H5Fclose(fileId);
  //    cerr <<"exiting fields::rstoreH5.\n";
  return(1);
}
#endif


int Fields::Restore(FILE *DMPFile)
{
  Scalar temp;
  int Jl,Kl;  // dimensions of the data in the dump file
  int jl,kl;  // indices over the data in the dump file
  
  
  XGRead(&Jl,4,1,DMPFile,"int");
  XGRead(&Kl,4,1,DMPFile,"int");
 
  if(!ElectrostaticFlag)
    {
      cout << "Restoring Fields:" << endl;
      cout << "Jl = " << Jl << endl;
      cout << "Kl = " << Kl << endl;

      // restoring fields
      for (jl=0; jl<=Jl;jl++)
        for (kl=0; kl<=Kl;kl++)
          {
            Scalar A1,A2;  // coordinates of the fields read in
            Vector2 Xloc;
            Vector3 E,B,EN,BN,Il;
            
            
            // read in the coordinates
            XGRead(&A1,ScalarInt,1,DMPFile,ScalarChar);
            XGRead(&A2,ScalarInt,1,DMPFile,ScalarChar);
	
            // find the grid locations here
            Xloc = grid->getNearestGridCoords(Vector2(A1,A2));
            
            
            //E
            XGVector3Read(E, DMPFile);
            //B
            XGVector3Read(B, DMPFile);
            /*		  //ENode
                          XGVector3Read(EN, DMPFile);  // not needed, recalculable.
                          //BNode
                          XGVector3Read(BN, DMPFile); */
            //J
            XGVector3Read(Il, DMPFile);
            if(Xloc.e1() >= 0) // valid location
              {
                int ji,ki;
                // the following operations make sure
                // stuff doesn't get rounded down when
                // converted from float to int
                ji = (int)(Xloc.e1()+0.5);
                ki = (int)(Xloc.e2()+0.5);
                
                intEdl[ji][ki]=E;
                intBdS[ji][ki]=B;
                /*ENode[ji][ki]=EN;
                  BNode[ji][ki]=BN; */ // Not needed: recalculable.
                I[ji][ki]=Il;
              }
          }
      XGRead(&temp, ScalarInt, 1, DMPFile, ScalarChar);
      if(temp>0){
        
        // red in the EMDamping stuff
        for (jl=0; jl<=Jl;jl++)
          for (kl=0; kl<=Kl;kl++)
            {
              Scalar A1,A2;  // coordinates of the fields read in
              Vector2 Xloc;
              Vector3 Emd;
              
              // read in the coordinates
              XGRead(&A1,ScalarInt,1,DMPFile,ScalarChar);
              XGRead(&A2,ScalarInt,1,DMPFile,ScalarChar);
        
	    
      
              // find the grid locations here
              Xloc = grid->getNearestGridCoords(Vector2(A1,A2));
	
              //Emd
              XGVector3Read(Emd, DMPFile);
              if(Xloc.e1() >= 0) // valid location
                {
                  int ji,ki;
                  // the following operations make sure
                  // stuff doesn't get rounded down when
                  // converted from float to int
                  ji = (int)(Xloc.e1()+0.5);
                  ki = (int)(Xloc.e2()+0.5);
		                    
                  intEdlBar[ji][ki]=Emd;
                }
            }			
        
        
        // if we've just turned EMDamping on.
        if (EMdamping != temp){
          cout << "**********************warning**********************" << endl;
          cout << "EMdamping value has been turned on, using damping." << endl;
          cout << "Weighted average is beening started with the instantaneous" << endl;
          cout << "value of the fields." << endl;
          for (int j=0; j<=J;j++)
            for (int k=0; k<=K;k++)
              intEdlBar[j][k] = intEdl[j][k];
        }
      }
    }
  
  if (BoltzmannFlag)
    theBoltzmann->Restore(DMPFile);
  
  return(1);
}


//
// Routines for moving fields around:  left, right, down, or up.
//

void Fields::sendFieldsAndShift(int dir, Boundary* bPtr){
//  int jp1; // temporary variable...
// Send cells over.  (JRC 22 Jan 00)
// This should send j=1/2 and j=1
  if(bPtr) bPtr->sendStripe(dir);

   //
   // Shift fields.  We do this first because the before-shift
   // edge cells were already sent over at the field solve.
   // I am not sure that this is correct.
   //
   switch(dir) {

      // Shift fields to the left
    case 1:
      {
        for(int j=0;j<J;j++) {
	  for(int k=0;k<=K;k++) {
	    ENode[j][k] = ENode[j+1][k];
	    BNode[j][k] = BNode[j+1][k];
	    BNodeDynamic[j][k] = BNodeDynamic[j+1][k];
	    BNodeStatic[j][k] = BNodeStatic[j+1][k];
	    I[j][k] = I[j+1][k];
	    intEdl[j][k] = intEdl[j+1][k];
	    intBdS[j][k] = intBdS[j+1][k];
            if(accumulator)
              accumulator[j][k] = accumulator[j+1][k];
            
            // stuff added by Bruhwiler (11/08/99) in attempt to prevent
            // moving window algorithm from crashing emdamping algorithm
	    intEdlBar[j][k] = intEdlBar[j+1][k];
	    intEdlPrime[j][k] = intEdlPrime[j+1][k];
	    intEdlBasePtr[j][k] = intEdlBasePtr[j+1][k];
	    Inode[j][k] = Inode[j+1][k];
	    rho[j][k] = rho[j+1][k];
	    backGroundRho[j][k] = backGroundRho[j+1][k];
	    DivDerror[j][k] = DivDerror[j+1][k];
	    Phi[j][k] = Phi[j+1][k];
	    PhiP[j][k] = PhiP[j+1][k];
	    Charge[j][k] = Charge[j+1][k];
	  }
        }
      }
      break;
      
      // Shift fields to the right
     case 2:
        {

           for(int j=J;j>1;j--) {
              for(int k=0;k<=K;k++) {
                 ENode[j][k] = ENode[j-1][k];
                 BNode[j][k] = ENode[j-1][k];
                 BNodeDynamic[j][k] = BNodeDynamic[j-1][k];
                 BNodeStatic[j][k] = BNodeStatic[j-1][k];
                 I[j][k] = I[j-1][k];
                 intEdl[j][k] = intEdl[j-1][k];
                 intBdS[j][k] = intBdS[j-1][k];
                 if(accumulator)
                   accumulator[j][k] = accumulator[j-1][k];
              }
           }
        }
       break;

       // Shift fields upward (smaller k)
     case 3:
        {
           for(int j=0;j<=J;j++) {
              for(int k=0;k<K;k++) {
                 ENode[j][k] = ENode[j][k+1];
                 BNode[j][k] = ENode[j][k+1];
                 BNodeDynamic[j][k] = BNodeDynamic[j][k+1];
                 BNodeStatic[j][k] = BNodeStatic[j][k+1];
                 I[j][k] = I[j][k+1];
                 intEdl[j][k] = intEdl[j][k+1];
                 intBdS[j][k] = intBdS[j][k+1];
                 if(accumulator)
                   accumulator[j][k] = accumulator[j][k+1];
              }
           }
        }
       break;

       // Shift fields downward (larger k)
     case 4:
        {
           for(int j=0;j<=J;j++) {
              for(int k=K;k>0;k--) {
                 ENode[j][k] = ENode[j][k-1];
                 BNode[j][k] = ENode[j][k-1];
                 BNodeDynamic[j][k] = BNodeDynamic[j][k-1];
                 BNodeStatic[j][k] = BNodeStatic[j][k-1];
                 I[j][k] = I[j][k-1];