fields.cpp

pic 模拟程序！面向对象

  for(j=0;j<J;j++)
     for(k=0;k<K;k++) {
        setIntEdl(j,k,1,
                  intEdl[j][k].e1() - (delPhi[j][k]-delPhi[j+1][k])   );
        setIntEdl(j,k,2,  
                  intEdl[j][k].e2() - (delPhi[j][k]-delPhi[j][k+1])   );
     }
}

void Fields::InodeToI(void)
{ int j,k;

  for(j=0;j<J;j++) for(k=0;k<K;k++)
     I[j][k] = 0.5*Vector3( 
                           ((j==0)?2:1)*Inode[j][k].e1() + ((j==J-1)?2:1)*Inode[j+1][k].e1(),
                           Inode[j][k].e2() + Inode[j][k+1].e2(),
                           2*Inode[j][k].e3());
}

//---------------------------------------------------------------------
//  Do a marder correction pass.  Actual algorithm is based on
//  Bruce Langdon's paper "On Enforcing Gauss' Law"

void Fields::MarderCorrection(Scalar dt)
{ int i,j,k;

  static Scalar** d=0;
  Vector2 **X=grid->getX();
  if(d==0) {

     d = new Scalar* [J+1];
     for(j=0;j<=J;j++)
        {
           d[j] = new Scalar[K+1];

           for(k=0;k<=K;k++) {
              if(j==J||k==K) d[j][k]=0; 
              else

                 d[j][k] = 0.5 * MarderParameter*iEPS0
                    /( 1.0/sqr( X[j+1][k].e1()-X[j][k].e1())
                      +1.0/sqr( X[j][k+1].e2()-X[j][k].e2()));
           }
        }
  }

  for(i=0;i<MarderIter;i++) {

     updateDivDerror();
     for(j=0;j<J;j++)
        for(k=0;k<K;k++) { 

           setIntEdl(j,k,1,intEdl[j][k].e1() - (d[j+1][k]*DivDerror[j+1][k] 
                                                - d[j][k]*DivDerror[j][k]));
           setIntEdl(j,k,2,intEdl[j][k].e2() - (d[j][k+1]*DivDerror[j][k+1] 
                                                - d[j][k]*DivDerror[j][k]));

        }
  }
}

void Fields::compute_iC_iL()
{
   Vector3** dl = grid->dl();
   Vector3** dlPrime = grid->dlPrime();
   Vector3** dS = grid->dS();
   Vector3** dSPrime = grid->dSPrime();
   Scalar epsjk, epsjmk, epsjkm, epsjmkm;

   for (int j=0; j<=J; j++)
      for (int k = 0; k<=K; k++){
         iL[j][k] = iMU0*dlPrime[j][k].jvDivide(dS[j][k]+Vector3(1e-30,1e-30,1e-30));
         epsjk = epsi[MIN(j,J-1)][MIN(k,K-1)];
         epsjmk = epsi[MAX(j-1, 0)][MIN(k,K-1)];
         epsjkm = epsi[MIN(j,J-1)][MAX(k-1, 0)];
         epsjmkm = epsi[MAX(j-1, 0)][MAX(k-1, 0)];
         iC[j][k] = dl[j][k].jvDivide(dSPrime[j][k].jvMult(
                                                           Vector3(.5*(epsjk+epsjkm),.5*(epsjk+epsjmk),.25*(epsjk+epsjmk+epsjkm+epsjmkm))));
      }
}

// morgan
void Fields::set_up_inverter(Grid* grid)
  throw(Oops){
   int nc1 = grid->getJ(); 
   int nc2 = grid->getK();

   // calculate uniform mesh widths for passing to Domain object
   Vector2 x;
   x = grid->getMKS(nc1,nc2);
   Scalar len1 = x.e1();
   Scalar len2 = x.e2();

   bool periodic1 = (grid->getPeriodicFlagX1() != 0); 
   try{
     dom = new Domain(grid,len1,len2,nc1,nc2,periodic1);
   }
   catch(Oops& oops2){
      oops2.prepend("Fields::set_up_inverter: Error: \n"); //Fields::initPhi
      throw oops2;
   }

   poisson = new Electrostatic_Operator(dom,epsi); 
   
} 



#ifdef HAVE_HDF5
int Fields::dumpH5(dumpHDF5 &dumpObj)
{
  
  // cerr << "Entered fields::dumpH5(dumpObj) .\n";
 
   hid_t fileId, groupId, datasetId,dataspaceId;
  herr_t status;
    hid_t scalarType = dumpObj.getHDF5ScalarType();
  hsize_t RANK  = 3;
//  hsize_t     dimsf[3];              /* dataset dimensions */
  
//  hssize_t start1[1]; /* Start of hyperslab */
//  hsize_t stride1[1]; /* Stride of hyperslab */
//  hsize_t count1[1];  /* Block count */
//  hsize_t block1[1];  /* Block sizes */

//  hssize_t start2[2]; /* Start of hyperslab */
//  hsize_t stride2[2]; /* Stride of hyperslab */
//  hsize_t count2[2];  /* Block count */
//  hsize_t block2[2];  /* Block sizes */
  
//  hssize_t start[3]; /* Start of hyperslab */
//  hsize_t stride[3]; /* Stride of hyperslab */
//  hsize_t count[3];  /* Block count */
//  hsize_t block[3];  /* Block sizes */


//       /* data to read */
//    int         i, j, num;
//    Scalar         val;

//    Scalar temp;
    Vector2 **X = grid->getX();
  //  Scalar x,y;


    fileId= H5Fopen(dumpObj.getDumpFileName().c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
    groupId = H5Gcreate(fileId,"Fields",0);
 
    if(!ElectrostaticFlag)
      {		  
	cout << "Dumping fields" << endl;
	cout << "J = " << J << endl;
	cout << "K = " << K << endl;
     
	hsize_t dim1 = 2;
	hsize_t dim3 = 3;
	//	hid_t posSlabId = H5Screate_simple(1, &dim1, NULL); // dimensions of position hyperslab
	//     * Create dataspace for the dataset in the file.
	hsize_t posDim[3];
	posDim[0] = J+1; posDim[1] = K+1; posDim[2] = 2;

	dataspaceId = H5Screate_simple(3, posDim, NULL);
	datasetId = H5Dcreate(groupId, "positions", scalarType, dataspaceId, H5P_DEFAULT);
	
	hsize_t fieldDim[3];
	fieldDim[0] = J+1; fieldDim[1] = K+1; fieldDim[2] = 3;
	//hid_t fieldSlabId = H5Screate_simple(1, &dim3, NULL);  // dimensions of field hyperslab

	hid_t EdlspaceId = H5Screate_simple(3, fieldDim, NULL);
	hid_t EdlsetId = H5Dcreate(groupId, "Edl", scalarType, EdlspaceId, H5P_DEFAULT);
	
	hid_t BdSspaceId = H5Screate_simple(3, fieldDim, NULL);
	hid_t BdSsetId = H5Dcreate(groupId, "BdS", scalarType, BdSspaceId, H5P_DEFAULT);

	hid_t IspaceId = H5Screate_simple(3, fieldDim, NULL);
	hid_t IsetId = H5Dcreate(groupId, "I", scalarType, IspaceId, H5P_DEFAULT);

	
	Scalar *posData = new Scalar[(J+1)*(K+1)*2]; //contiguous values needed
	

	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++)
	    {
	     //   if(j*(K+1) +k  < 10){
//  		cerr << X[j][k].e1()<<"\n";
//  		cerr << X[j][k].e2()<<"\n";
//  	      }

	      posData[(K+1)*2*j + k*2]= X[j][k].e1();
	      posData[(K+1)*2*j + k*2 + 1]= X[j][k].e2();
	    }


	status = H5Dwrite(datasetId, scalarType, H5S_ALL,H5S_ALL , H5P_DEFAULT, posData);
	delete[] posData;

	Scalar *data = new Scalar[(J+1)*(K+1)*3];
    
	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++)
	    {
	      data[(K+1)*3*j + k*3] = intEdl[j][k].e1();
	      data[(K+1)*3*j + k*3 + 1] = intEdl[j][k].e2();
	      data[(K+1)*3*j + k*3 + 2] = intEdl[j][k].e3();
	    }

	status = H5Dwrite(EdlsetId, scalarType, H5S_ALL, EdlspaceId, H5P_DEFAULT, data);

	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++)
	    {
	      data[(K+1)*3*j + k*3]= intBdS[j][k].e1();
	      data[(K+1)*3*j + k*3 + 1]= intBdS[j][k].e2();
	      data[(K+1)*3*j + k*3 + 2]= intBdS[j][k].e3();

	      


   	    }
	status = H5Dwrite(BdSsetId, scalarType, H5S_ALL, BdSspaceId, H5P_DEFAULT, data);

	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++)
	    {
	      data[(K+1)*3*j + k*3]= I[j][k].e1();
	      data[(K+1)*3*j + k*3 + 1]= I[j][k].e2();
	      data[(K+1)*3*j + k*3 + 2]= I[j][k].e3();
	    }
	status = H5Dwrite(IsetId, scalarType, H5S_ALL, IspaceId, H5P_DEFAULT, data);
       
	
	
	status = H5Sclose(dataspaceId); 
	status = H5Sclose(EdlspaceId); 
	status = H5Sclose(BdSspaceId); 
	status = H5Sclose(IspaceId); 
	/*
	 * Close datasets.
	 */
	status = H5Dclose(datasetId);
	status = H5Dclose(EdlsetId);
	status = H5Dclose(BdSsetId);
	status = H5Dclose(IsetId);
	

    // dump EM damping stuff
	
// write EMdamping flag as attribute

	hsize_t     dims;
//	int         attrData;
	hid_t  attrdataspaceId,attributeId;
	dims = 1;
	attrdataspaceId = H5Screate_simple(1, &dims, NULL);
	attributeId = H5Acreate(groupId, "EMdampingFlag",scalarType, attrdataspaceId, H5P_DEFAULT);
  
	status = H5Awrite(attributeId, scalarType, &EMdamping);
	status = H5Aclose(attributeId);
	status = H5Sclose(attrdataspaceId);



        EdlspaceId = H5Screate_simple(3, fieldDim, NULL);
	EdlsetId = H5Dcreate(groupId, "EdlBar", scalarType, EdlspaceId, H5P_DEFAULT);

	

	
	if((EMdamping)>0)
	  for (int j=0; j<=J;j++)
	    for (int k=0; k<=K;k++) {

	      data[(K+1)*3*j + k*3] = intEdlBar[j][k].e1();
	      data[(K+1)*3*j + k*3 + 1] = intEdlBar[j][k].e2();
	      data[(K+1)*3*j + k*3 + 2] = intEdlBar[j][k].e3();
	    }
	
	status = H5Dwrite(EdlsetId, scalarType, H5S_ALL, EdlspaceId, H5P_DEFAULT, data);
	status = H5Sclose(EdlspaceId); 
	status = H5Dclose(EdlsetId);
      }
    /*
     * Close the file.
     */
    status = H5Gclose(groupId);
    status = H5Fclose(fileId);
    // dump boltzmann

   if (BoltzmannFlag)
     theBoltzmann->dumpH5(dumpObj);
    
	return(1);
}

#endif

int Fields::Dump(FILE *DMPFile)
{
//   Scalar temp;
   Vector2 **X = grid->getX();
   Scalar x,y;
//   char theDumpFileName[256];
   char *fdname ={"fdmp"};

   XGWrite(&J,4,1,DMPFile,"int");
   XGWrite(&K,4,1,DMPFile,"int");
   
   if(!ElectrostaticFlag)
      {		  
      // ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.
/*
       #ifdef MPI_VERSION
           extern int MPI_RANK;
           sprintf(theDumpFileName,"%s.%d",fdname,MPI_RANK);
       #else
           sprintf(theDumpFileName,"%s",fdname);
       #endif
       ofstream dataField(theDumpFileName, ios::app);

*/
        cout << "Dumping fields" << endl;
        cout << "J = " << J << endl;
        cout << "K = " << K << endl;
#ifdef MPI_VERSION
      extern int MPI_RANK;
      cout << "MPI rank " << MPI_RANK << endl;
#endif

      cout << "\nField positions:\n";
      
      for (int j=0; j<=J;j++)
	    for (int k=0; k<=K;k++)
	       {
		  x = X[j][k].e1();
		  y = X[j][k].e2();
		  XGWrite(&x,ScalarInt,1,DMPFile,ScalarChar);
		  XGWrite(&y,ScalarInt,1,DMPFile,ScalarChar);

		  //E
		  XGVector3Write(intEdl[j][k], DMPFile);
		  //B
		  XGVector3Write(intBdS[j][k], DMPFile);
		  //Enode
		  /*XGVector3Write(ENode[j][k], DMPFile); // not needed, recalculable.
                    //Bnode
                    XGVector3Write(BNode[j][k], DMPFile); */		  
		  //J
		  XGVector3Write(I[j][k], DMPFile);
/*
	          //TEXT
		  dataField << x <<"  "<< y <<"  "<< ENode[j][k].e1() <<"  "<<
                           ENode[j][k].e2()<<"  "<<ENode[j][k].e3()<<endl;
*/
	       }
//        dataField.close();
	 XGWrite(&EMdamping, ScalarInt, 1, DMPFile, ScalarChar);
	 if((EMdamping)>0)
	    for (int j=0; j<=J;j++)
	       for (int k=0; k<=K;k++) {
		  x = X[j][k].e1();
		  y = X[j][k].e2();
		
		  XGWrite(&x,ScalarInt,1,DMPFile,ScalarChar);
		  XGWrite(&y,ScalarInt,1,DMPFile,ScalarChar);

		  XGVector3Write(intEdlBar[j][k], DMPFile);
	       }
      }
   if (BoltzmannFlag)
      theBoltzmann->Dump(DMPFile);

// ~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.~.

   return (1);

}

int Fields::Restore_2_00(FILE *DMPFile)
{
   Scalar temp;

   if(!ElectrostaticFlag)
      {
         // restoring fields
         for (int j=0; j<=J;j++)
            for (int k=0; k<=K;k++)
               {
                  //E
                  XGVector3Read(intEdl[j][k], DMPFile);
                  //B
                  XGVector3Read(intBdS[j][k], DMPFile);
                  //J
                  XGVector3Read(I[j][k], DMPFile);

               }
         XGRead(&temp, ScalarInt, 1, DMPFile, ScalarChar);
         if(temp>0){
            if(EMdamping==temp)
               for (int j=0; j<=J;j++)
                  for (int k=0; k<=K;k++)
                     XGVector3Read(intEdlBar[j][k], DMPFile);
            else if(EMdamping>=0)
               {
                  cout << "**********************warning**********************" << endl;
                  cout << "EMdamping value has been changed, using new value." << endl;
                  cout << "Keeping old weighted average." << endl;
                  for (int j=0; j<J;j++)
                     for (int k=0; k<K;k++)
                        XGVector3Read(intEdlBar[j][k], DMPFile);
               }
         }
         else 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)