ptclgrp.cpp

pic 模拟程序！面向对象

 
  int vary_np2c = 0;
  int speciesID = get_speciesID();
//  Scalar temp;
  Vector2 xMKS;
//  char theDumpFileName[256];
  char *fdname ={"pdmp"};

  hid_t fileId, groupId, datasetId,dataspaceId;
  herr_t status;
  hsize_t rank; 
  hsize_t     dimsf[2];   
  hsize_t     dims;
//  int         attrData;
  hid_t  attrdataspaceId,attributeId;
#ifdef NEW_H5S_SELECT_HYPERSLAB_IFC
  hsize_t start[2]; /* Start of hyperslab */
#else
  hssize_t start[2]; /* Start of hyperslab */
#endif
  hsize_t count[2];  /* Block count */

  Scalar (*qData) = new Scalar[n];
  hid_t scalarType = dumpObj.getHDF5ScalarType();
 
  if (qArray)	vary_np2c = 1;
 
 // Open the hdf5 file with write access
  fileId = H5Fopen(dumpObj.getDumpFileName().c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
  groupId = H5Gcreate(fileId,groupName.c_str(),0);

 
  // dump speciesID as attribute
   dims = 1;
   attrdataspaceId = H5Screate_simple(1, &dims, NULL);
   attributeId = H5Acreate(groupId, "speciesID",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
   status = H5Awrite(attributeId,H5T_NATIVE_INT , &speciesID);
   status = H5Aclose(attributeId);
   status = H5Sclose(attrdataspaceId);
   
 // dump vary_np2c as attribute
   dims = 1;
   attrdataspaceId = H5Screate_simple(1, &dims, NULL);
   attributeId = H5Acreate(groupId, "vary_np2c",H5T_NATIVE_INT, attrdataspaceId, H5P_DEFAULT);
   status = H5Awrite(attributeId,H5T_NATIVE_INT , &vary_np2c);
   status = H5Aclose(attributeId);
   status = H5Sclose(attrdataspaceId);

 // dump np2c as attribute
   dims = 1;
   attrdataspaceId = H5Screate_simple(1, &dims, NULL);
   attributeId = H5Acreate(groupId, "np2c",scalarType, attrdataspaceId, H5P_DEFAULT);
   status = H5Awrite(attributeId, scalarType, &np2c0);
   status = H5Aclose(attributeId);
   status = H5Sclose(attrdataspaceId);
 
   rank = 2;
   dimsf[0] = n;
   dimsf[1] = 5;
   dataspaceId = H5Screate_simple(rank, dimsf, NULL); 
   string datasetName = "ptcls";  
   datasetId = H5Dcreate(groupId, datasetName.c_str(), scalarType, dataspaceId,
			H5P_DEFAULT);


   hsize_t dim5 = 5;
   hid_t posSlabId = H5Screate_simple(1, &dim5, NULL);

   for(int i=0;i<n;i++){		
     
     xMKS = grid->getMKS(x[i]);

     Scalar pos[5] = {xMKS.e1(),xMKS.e2(),u[i].e1(),u[i].e2(),u[i].e3()};  
     start[0]  = 0;  count[0]  = 5;
     status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, start, NULL, count, NULL);
     start[0]  = i; start[1]  = 0;
     count[0]  = 1; count[1]  = 5;     
     status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, start, NULL, count, NULL);
     status = H5Dwrite(datasetId, scalarType, posSlabId, dataspaceId, H5P_DEFAULT, pos);
     
   }


   H5Sclose(posSlabId);
   H5Sclose(dataspaceId);
   H5Dclose(datasetId);


if (qArray)
  {
// write charge   
    rank = 1;
    dims = n;
    dataspaceId = H5Screate_simple(rank, &dims, NULL); 
    string datasetName = "Q";  
    datasetId = H5Dcreate(groupId, datasetName.c_str(), scalarType, dataspaceId,
			  H5P_DEFAULT);

    status = H5Dwrite(datasetId, scalarType, H5S_ALL, H5S_ALL,
		      H5P_DEFAULT, qArray); 

    H5Sclose(dataspaceId);
    H5Dclose(datasetId);
  }
 
 H5Gclose(groupId);
 H5Fclose(fileId);
 return (1);
  
}
#endif

int ParticleGroup::Dump(FILE *DMPFile,Grid *grid) 
{
	int vary_np2c = 0;
	int speciesID = get_speciesID();
//	Scalar temp;
	Vector2 xMKS;
//        char theDumpFileName[256];
        char *fdname ={"pdmp"};
/*
//reg
        #ifdef MPI_VERSION
	   extern int MPI_RANK;
           sprintf(theDumpFileName,"%s.%d",fdname,MPI_RANK);
	#else
           sprintf(theDumpFileName,"%s",fdname);
        #endif
//
	ofstream dataPart(theDumpFileName, ios::app);
*/

	XGWrite(&speciesID, 4, 1, DMPFile, "int");
	if (qArray)	vary_np2c = 1;
	XGWrite(&vary_np2c, 4, 1, DMPFile, "int");
	XGWrite(&np2c0, ScalarInt, 1, DMPFile, ScalarChar);
	XGWrite(&n, 4, 1, DMPFile, "int");

	if (qArray)
	   {
		for(int i=0;i<n;i++){
			XGWrite(&qArray[i], ScalarInt, 1, DMPFile, ScalarChar);
			xMKS = grid->getMKS(x[i]);
			XGVector2Write(xMKS, DMPFile);
			XGVector3Write(u[i], DMPFile);

		       
/*
			dataPart << xMKS.e1()<<"  " << xMKS.e2() 
				 <<"   "<< u[i].e1() << endl;
*/
		}
	   }
	else{
		for(int i=0;i<n;i++){
			xMKS = grid->getMKS(x[i]);
			XGVector2Write(xMKS, DMPFile);
			XGVector3Write(u[i], DMPFile);
			//			cout << xMKS.e1() << " "<<xMKS.e2()<<endl;

/*
			dataPart << xMKS.e1()<<"  " << xMKS.e2() 
				 <<"   "<< u[i].e1() << endl;
*/

		}
	}
//	dataPart.close();

	return (1);
} 


#ifdef HAVE_HDF5
int ParticleGroup::restoreH5(dumpHDF5 &restoreObj, Grid *grid, Scalar np2cFactor,string whichGroup)
{
  // cerr << "entering ParticleGroup::restoreH5()\n";

  Scalar xmin, xmax,ymin,ymax;
  Vector2 xin;
  Vector3 uin;
  Scalar qin;
  int i;
  int incount;
  int realcount=0;
  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,attrId;
#ifdef NEW_H5S_SELECT_HYPERSLAB_IFC
  hsize_t offset[2]; /* Offset of hyperslab */
#else
  hssize_t offset[2]; /* Offset of hyperslab */
#endif
  hsize_t count[2];  /* Block count */


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

  xmin = grid->getX()[0][0].e1();
  ymin = grid->getX()[0][0].e2();
  xmax = grid->getX()[grid->getJ()][grid->getK()].e1();
  ymax = grid->getX()[grid->getJ()][grid->getK()].e2();

  // read Q values
 
  // read position values
  datasetId = H5Dopen(groupId,"ptcls");
  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);
  incount = sizes[0];
  hsize_t dim = sizes[1];  // should be 5
  hid_t posSlabId =  H5Screate_simple(1, &dim, NULL);
  
  Scalar *ptclData = new Scalar[dim];


  if(qArray){
    hid_t QsetId = H5Dopen(groupId,"Q");
    hid_t QspaceId = H5Dget_space(QsetId );


    for(i=0;i<incount;i++){
      offset[0]  = 0;  count[0]  = dim;
      status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, offset, NULL, count, NULL);
      offset[0]  = i; offset[1]  = 0;
      count[0]  = 1; count[1] = dim;    
      status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, offset, NULL, count, NULL);
      status = H5Dread(datasetId, scalarType, posSlabId, dataspaceId, H5P_DEFAULT, ptclData);
      
      xin.set_e1(ptclData[0]);
      xin.set_e2(ptclData[1]);
      
      //  cout << ptclData[2] <<" " <<ptclData[3]<<" "<<ptclData[4]<<"\n";
      
      
      uin.set_e1(ptclData[2]);
      uin.set_e2(ptclData[3]);
      uin.set_e3(ptclData[4]);
      
      offset[0]  = 0;  count[0]  = 1;
      status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, offset, NULL, count, NULL);
      offset[0]  = i; offset[1]  = 0;
      count[0]  = 1; count[1] = 1;    
      status = H5Sselect_hyperslab(QspaceId, H5S_SELECT_SET, offset, NULL, count, NULL);
      status = H5Dread(QsetId, scalarType, posSlabId, QspaceId, H5P_DEFAULT, &qin);
      
      // check to see if inbounds
      if(xin.e1() > xmin && xin.e1() < xmax &&
	 xin.e2() > ymin && xin.e2() < ymax) {
	qArray[realcount]=qin*np2cFactor;
	x[realcount] = grid->getGridCoords(xin);
	u[realcount] = uin;
	realcount++;
      }
    }
    
    status =    H5Sclose(QspaceId);
    status =    H5Dclose(QsetId);
    
  }
  else {
    
    //      cout <<"particles at restore:\n";
    //   cerr << "incount"<<incount<<"\n";
    
    for(i=0;i<incount;i++){
      
      offset[0]  = 0;  count[0]  = dim;
      status = H5Sselect_hyperslab(posSlabId, H5S_SELECT_SET, offset, NULL, count, NULL);
      offset[0]  = i; offset[1]  = 0;
      count[0]  = 1; count[1] = dim;    
      status = H5Sselect_hyperslab(dataspaceId, H5S_SELECT_SET, offset, NULL, count, NULL);
      status = H5Dread(datasetId, scalarType, posSlabId, dataspaceId, H5P_DEFAULT, ptclData);
      
      xin.set_e1(ptclData[0]);
      xin.set_e2(ptclData[1]);
      
      uin.set_e1(ptclData[2]);
      uin.set_e2(ptclData[3]);
      uin.set_e3(ptclData[4]);
      
      // check to see if inbounds
      if(xin.e1() > xmin && xin.e1() < xmax &&
	 xin.e2() > ymin && xin.e2() < ymax) {
	x[realcount] = grid->getGridCoords(xin);
	u[realcount] = uin;
	realcount++;
      }
    }
    
  }
  n = realcount;
 
  status =    H5Sclose(dataspaceId);
  status =    H5Dclose(datasetId);
 
  status = H5Sclose(posSlabId);
  H5Gclose(groupId);
  H5Fclose(fileId);

  delete [] sizes;
  return (1);
}
#endif

int ParticleGroup::Restore(FILE *DMPFile,Grid *grid, Scalar np2cFactor) 
{
	Scalar xmin, xmax,ymin,ymax;
	int i;
	int incount;
	int realcount=0;
	Vector2 xin;
	Vector3 uin;
	Scalar qin;

	xmin = grid->getX()[0][0].e1();
	ymin = grid->getX()[0][0].e2();
	xmax = grid->getX()[grid->getJ()][grid->getK()].e1();
	ymax = grid->getX()[grid->getJ()][grid->getK()].e2();

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

	if (qArray)
	  for(i=0;i<incount;i++){
			XGRead(&qin, ScalarInt, 1, DMPFile, ScalarChar);
			XGVector2Read(xin, DMPFile);
			XGVector3Read(uin, DMPFile);

			// check to see if inbounds
			if(xin.e1() > xmin && xin.e1() < xmax &&
				xin.e2() > ymin && xin.e2() < ymax) {
				qArray[realcount]=qin*np2cFactor;
				x[realcount] = grid->getGridCoords(xin);
				u[realcount] = uin;
				realcount++;
			}
		}
	else
		for(i=0;i<incount;i++){
			XGVector2Read(xin, DMPFile);
			XGVector3Read(uin, DMPFile);

			// check to see if inbounds
			if(xin.e1() > xmin && xin.e1() < xmax &&
				xin.e2() > ymin && xin.e2() < ymax) {
				x[realcount] = grid->getGridCoords(xin);
				u[realcount] = uin;
				realcount++;
			}
		}

	n = realcount;
	return(1);
}

int ParticleGroup::Restore_2_00(FILE *DMPFile) 
{

	XGRead(&n, 4, 1, DMPFile, "int");
	if (qArray)
		for(int i=0;i<n;i++){
			XGRead(&qArray[i], ScalarInt, 1, DMPFile, ScalarChar);
			XGVector2Read(x[i], DMPFile);
			XGVector3Read(u[i], DMPFile);
		}
	else
		for(int i=0;i<n;i++){
			XGVector2Read(x[i], DMPFile);
			XGVector3Read(u[i], DMPFile);
		}

	return(1);

}


void ParticleGroup::shiftParticles(Fields& fields, Vector2 shift) {

  Vector2*	x = ParticleGroup::x;
  Vector3*	u = ParticleGroup::u;
  Grid *grid = fields.get_grid();
  Boundary *bPtr;
  Scalar smallTest = 1.0e-06;
  Scalar x_e1_tmp;

  for(int i=0;i<n; i++,x++,u++) {

// Shift the particle
    *x += shift;

// Start with no boundary
    bPtr = 0;
    bool removeParticle = false;

// Take care of particles near the boundary.
// Unfortunately, this kind of nonsense is required to prevent
//    particles from being caught "on" the edge of a boundary.
// First, consider particles that should be left behind.
    x_e1_tmp = x->e1();
    if( ( x_e1_tmp >= 0. )  &&  ( x_e1_tmp < smallTest ) ) {
      x->set_e1(smallTest);
		}
// Next, consider particles that need to be shifted.
    else if( ( x_e1_tmp <= 0. )  &&  ( x_e1_tmp > -smallTest ) ) {
      x->set_e1(-smallTest);
		}

// Now do the shift, assuming a rightward moving window.
    if( (shift.e1() == -1.) && ( x->e1() < 0. ) ) {
// Peter's original code:
//      if( ( bPtr=grid->GetBC2()[(int)(*x).e1()][(int)(*x).e2()] ) == 0 )
//                bPtr = grid->GetBC1()[(int)(*x).e1()][(int)(*x).e2()];
// relied on the fact that int of a negative number gives zero.
      bPtr = grid->GetBC2()[0][(int)(*x).e2()];
      if(!bPtr) {
        static bool printed = false;
        if(!printed) cout << "No boundary in shiftParticles!" << endl;
        printed = true;
      }
      removeParticle = true;
    }

// Give the particle to the boundary
    if(bPtr) bPtr->collectShiftedParticle( new Particle(*x, *u, species, 
		get_np2c(i), (BOOL)(qArray!=0) ) );
    if(removeParticle){

// Decrement number of particles
      n--;
      if (i == n) break;          // last particle in array?
      *x = ParticleGroup::x[n];   // move last particle into open slot
      *u = ParticleGroup::u[n];
      if (qArray) qArray[i] = qArray[n];

// Decrement counters as new particle moved here
      i--;
      x--;
      u--;

    }

  }

}