diagn.cpp

pic 模拟程序！面向对象

				for(thisSpc->ngroups=0,pgscan.restart();!pgscan.Done();pgscan++,thisSpc->ngroups++){
					pgnumparticles=(int)pgscan.current()->get_n();
					thisSpc->nparticles+=pgnumparticles;
				}

				// Skip particles in the physics so we don't overflow the
				// particle array in the diagnostics
				skip = thisSpc->nparticles/thisSpc->memory_size +1;
				thisSpc->nparticles_plot = thisSpc->nparticles/skip;
#ifdef MPI_VERSION
				Scalar l_nparticles = thisSpc->nparticles;
				Scalar t_nparticles = 0;
				MPI_Reduce ( &l_nparticles, &t_nparticles,1,MPI_Scalar,MPI_SUM,0,XOOPIC_COMM);
				if(MPI_RANK==0)
				  thisSpc->nparticles = (int)t_nparticles;
#endif /* MPI_VERSION */
				
				for(pgscan.restart();!pgscan.Done();pgscan++) {
					int npart =(int) pgscan.current()->get_n();
					for(j=0;j<npart&&i<NPARTICLES;j+=skip,i++) {
							Vector2 *xarr=pgscan.current()->get_x();
							Vector3 *varr=pgscan.current()->get_u();
							Vector2 x = theSpace->getMKS(xarr[j]);
							Vector3 v = varr[j];
							
							thisSpc->x1positions[i]=x.e1();
							thisSpc->x2positions[i]=x.e2();
							thisSpc->x1velocities[i] = v.e1();
							thisSpc->x2velocities[i] = v.e2();
							thisSpc->x3velocities[i] = v.e3();
						}
				}
			}
		}
		
	}
#endif
}

void Diagnostics::UpdatePostDiagnostics() {
  int j,k;
  int jm = theSpace->getJ();
  int km = theSpace->getK();

	if (boltzmannFlag)
		theSpace->CopyBoltzmannRho();

  //  initialize time history varaibles to zero
  Uetot=0; Ubtot=0;
  if(AllDiagnostics) {  //collect diagnostics at all?
		Grid* grid = theSpace->get_grid();
		theSpace->updateDivDerror();

		if(EnerPlots){
			for( j=0;j<=jm;j++)
				for( k=0;k<=km;k++){

					/* get the E energy density */
                                  Ue[j][k] = 0.5*(theSpace->get_fields()->get_epsi(MIN(j,jm-1),MIN(k,km-1))+
                                                  theSpace->get_fields()->get_epsi(MAX(j-1,0),MIN(k,km-1))+
                                                  theSpace->get_fields()->get_epsi(MIN(j,jm-1),MAX(k-1,0))+
                                                  theSpace->get_fields()->get_epsi(MAX(j-1,0),MAX(k-1,0)))/4.0
                                    *(E[j][k]*E[j][k]);
					/* get the B energy density */
					Ub[j][k]=sqr(BDynamic[j][k].e1());
					Ub[j][k]+=sqr(BDynamic[j][k].e2());
					Ub[j][k]+=sqr(BDynamic[j][k].e3());
					Ub[j][k]*=0.5*iMU0;

#ifdef DEBUG_PHI
					if(phi[j][k] > 1e-20)
						//absolute error (uncomment one)
					phi_err[j][k] = fabs((phi[j][k] - (sin( M_PI * grid->getX()[j][k].e1()) 
						* sin(M_PI*grid->getX()[j][k].e2()))));

						//relative error
//					phi_err[j][k] = fabs((phi[j][k] - (sin( M_PI * grid->getX()[j][k].e1()) 
//						* sin(M_PI*grid->getX()[j][k].e2())))/phi[j][k]);
					else phi_err[j][k] = 0;
#endif
                                          //collect totals of energy for volumetric data
					Uetot+=Ue[j][k]*CellVolumes[j][k];
					Ubtot+=Ub[j][k]*CellVolumes[j][k];
				}
			
			// This section makes the poyting vector on the interior points. 
			// To get the edge points it is a little harder.
			if (!electrostaticFlag){
				for( j=0;j<=(jm-1);j++)
					for( k=1;k<=(km-1);k++)
						S_array[j][k] = ((E[j][k].cross(BDynamic[j][k]))).
							jvMult(grid->dS(j,k))*iMU0;
			}
		}

		history();
		
  }
} 
#ifdef BENCHMARK
//  This function writes a 'trace', 
//  the first component of the electric field
//  so that you can check the output of a simulation
//  which has no GUI.

void write_validation() {
  FILE *trace_file;
  int J = theSpace->getJ();
  int K = theSpace->getK();

  if((trace_file=fopen("trace.dat","w"))==NULL) {
    stringstream ss (stringstream::in | stringstream::out);
    ss<< "diagn::write_validation: Error: \n"<<
      "can not open trace.dat \n"<<endl;
    
    std::string msg;
    ss >> msg;
    Oops oops(msg);
    throw oops;    // exit()  not called

  }
  for(int j=0;j<J;j++) 
    for(int k=0;k<K;k++) 
			fprintf(trace_file,"%10.4g\n",E[j][k].e1());

  fclose(trace_file);
}
#endif


#ifdef HAVE_HDF5
void Diagnostics::dumpGridH5(dumpHDF5 *H5DumpObj)  
{
  
  //cerr << "Dumping grid information to diagnostic dump file.\n";
  int rank = 2;
  int *size;
  int J,K;
  
// dump grid information
  J = theSpace->getJ();
  K = theSpace->getK();
  rank = 1;
  size = new int[rank];
  
  size[0] = J+1;
  H5DumpObj->writeSimple("x1array",x1_array,rank,size);
  
  size[0] = K+1;
  H5DumpObj->writeSimple("x2array",x2_array,rank,size);
  
  delete[] size;
  return;
}

void Diagnostics::dumpAllDiagsH5(oopicList<dumpHDF5> *H5DiagDumpObjs) 
{

// create iterator over dump objects, increment at end of loop
  oopicListIter<dumpHDF5> nObj(*H5DiagDumpObjs);
// for now assume all diagnostics output to same file, so only dump grid info once
  nObj.restart();
  dumpGridH5(nObj.current() );
 
    for (nObj.restart(); !nObj.Done(); nObj++){
//      cerr << "datasetName "<<(nObj.current())->datasetName<<"\n";
      dumpDiagH5(nObj.current());
    }

  return;
}


void Diagnostics::dumpDiagH5(dumpHDF5 *H5DumpObj) 
{

//  hid_t fileId, groupId, datasetId,dataspaceId;
//  herr_t status;
      
  string diagName = H5DumpObj->datasetName;
  //  cerr << "Dumping "<<diagName<<" diagnostic.\n";

    Scalar *data; 
    int rank = 2;
    int *size;
    int J,K;

    //  dumpGridH5(H5DumpObj);

    J = theSpace->getJ();
    K = theSpace->getK();
  
    if(diagName=="E1" || diagName=="E2" || diagName=="E3" || 
       diagName=="B1" || diagName=="B2" || diagName=="B3" || 
       diagName=="I1" || diagName=="I2" || diagName=="I3" ){
      rank = 3;
      size = new int[rank];
      
      size[0] = J+1; size[1] = K+1; size[2] = 1;
      data = new Scalar[(J+1)*(K+1)*1];
      
    
      
      
      if(diagName=="E1"){
	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++){
	    data[j*(K+1) + k] = E[j][k].e1();
	  }
      }else
	if(diagName=="E2"){
	  for (int j=0; j<=J;j++)
	    for (int k=0; k<=K;k++){
	      data[j*(K+1) + k] = E[j][k].e2();
	    }
	  
      }else  
      if(diagName=="E3"){
	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++){
	    data[j*(K+1) + k] = E[j][k].e3();
	  }

      }else
      if(diagName=="B1"){
	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++){
	    data[j*(K+1) + k] = B[j][k].e1();
	  }

      }else
	if(diagName=="B2"){
	  for (int j=0; j<=J;j++)
	    for (int k=0; k<=K;k++){
	    data[j*(K+1) + k] = B[j][k].e2();
	    }
	  
	}else
	  if(diagName=="B3"){
	    for (int j=0; j<=J;j++)
	      for (int k=0; k<=K;k++){
		data[j*(K+1) + k] = B[j][k].e3();
	      }
	  }else
      if(diagName=="I1"){
	for (int j=0; j<=J;j++)
	  for (int k=0; k<=K;k++){
	    data[j*(K+1) + k] = I[j][k].e1();
	  }

      }else
	if(diagName=="I2"){
	  for (int j=0; j<=J;j++)
	    for (int k=0; k<=K;k++){
	    data[j*(K+1) + k] = I[j][k].e2();
	    }
	  
	}else
	  if(diagName=="I3"){
	    for (int j=0; j<=J;j++)
	      for (int k=0; k<=K;k++){
		data[j*(K+1) + k] = I[j][k].e3();
	      }
	  }    
    
      // H5DumpObj->createFile();
      H5DumpObj->writeSimple(diagName,data,rank,size);
      
      delete[] size;
      delete[] data;
      return;
   
    }else
      if(diagName=="NGD"){
	rank = 3;
	size = new int[rank];
	
	size[0] = J; size[1] = K; size[2] = 1;
	data = new Scalar[(size[0])*(size[1])*(size[2])];
	
	
	if(diagName=="NGD"){
	  
	  // for now take the first NDG data pointer
	  oopicListIter<NGD> NGDIter(*ptrNGDList);
	  NGDIter.restart();  
	  
	  for (int j=0; j<J;j++)
	    for (int k=0; k<K;k++){
	      data[j*(K) + k] = ((NGDIter.current())->getNGDdata())[j][k];
	    }
	}
	  
	  
	H5DumpObj->writeSimple(diagName,data,rank,size);
	
	  
	// write the grid points for NGD
	string datasetName;
	datasetName ="/NGDx1array";
	rank = 1;
	size[0] = J;
	H5DumpObj->writeSimple(datasetName,x1_arrayNGD,rank,size);
	  
	datasetName = "/NGDx2array";
	size[0] = K;
	H5DumpObj->writeSimple(datasetName,x2_arrayNGD,rank,size);	
	
	
	
	delete[] size;
	delete[] data;
	return;
      }else
	for(int isp=0;isp<number_of_species;isp++) {
// loop over all species, dumping number density for each
// if it exists 
	  cerr << "diagName : "<<diagName<<"\n";
	  string nm = theSpecies[isp].name;
	  nm = nm + "_numberDensity";	
	  if(diagName == nm){
	    rank = 3;
	    size = new int[rank];
	    size[0] = J+1; size[1] = K+1; size[2] = 1;
	    data = new Scalar[(J+1)*(K+1)*1];
	    for (int j=0; j<=J;j++)
	      for (int k=0; k<=K;k++){
		data[j*(K+1) + k] = rho_species[isp][j][k];
	      }

	    cerr << "data: "<< data[size[0]-1]<<"\n";

	    H5DumpObj->writeSimple(nm,data,rank,size);
	   
	  
	  delete[] size;
	  delete[] data;
	  return;
	  }
	}
    
    
    
    DiagList *dlist = theSpace->getDiagList();
    oopicListIter<Diag> nextd(*dlist);

    int iii =0;
    for(nextd.restart();!nextd.Done(); nextd++){ 
      
      cerr << "names of diag "<<iii<<": "<< nextd.current()->getVarName().str()<<"\n";
      iii++;
}

    for(nextd.restart();!nextd.Done(); nextd++){ 
      cerr << "next names of diag "<<iii<<": "<< nextd.current()->getVarName().str()<<"\n";
      if( nextd.current()->getVarName().str()==diagName){
	nextd.current()->dumpDiagH5(H5DumpObj);
      }
    }
    return;
    
    
    cerr << "Input file calls for dump of " << diagName
	 << " diagnostic, which does not exist.\n\n";
    return;
}



void Diagnostics::dumpH5(dumpHDF5 &dumpObj) 
{
  //   cerr << "Entered Diagnostics::dumpH5(dumpObj) .\n";
 
  //   hid_t fileId, groupId, datasetId,dataspaceId;
//     herr_t status;
//     hid_t scalarType = dumpObj.getHDF5ScalarType();


//     fileId= H5Fopen(dumpObj.getDumpFileName().c_str(), H5F_ACC_RDWR, H5P_DEFAULT);
//     groupId = H5Gcreate(fileId,"Diagnostics",0);
//     int J = theSpace->getJ();
//     int K = theSpace->getK();
//     hsize_t fieldDim[3];
//     fieldDim[0] = J+1; fieldDim[1] = K+1; fieldDim[2] = 1;
   
  
//     hid_t EspaceId = H5Screate_simple(3, fieldDim, NULL);
//     hid_t EsetId = H5Dcreate(groupId, "E1", scalarType, EspaceId, H5P_DEFAULT);

//     Scalar *Edata = new Scalar[(J+1)*(K+1)*1];
   
//     for (int j=0; j<=J;j++)
//       for (int k=0; k<=K;k++){
//         Edata[j*(K+1) + k] = E[j][k].e1();
//         //   Edata[j*(K+1)*3 + k*3 + 1] =  E[j][k].e2();
//         //  Edata[j*(K+1)*3 + k*3 + 2] = E[j][k].e3();
//       }
//     status = H5Dwrite(EsetId, scalarType, H5S_ALL, EspaceId, H5P_DEFAULT, Edata);
//     status = H5Dclose(EsetId);
   
//     EsetId = H5Dcreate(groupId, "E2", scalarType, EspaceId, H5P_DEFAULT);
//     for (int j=0; j<=J;j++)
//       for (int k=0; k<=K;k++){
//         Edata[j*(K+1) + k] =  E[j][k].e2();
//        //  Edata[j*(K+1)*3 + k*3 + 2] = E[j][k].e3();
//  	    }
//     status = H5Dwrite(EsetId, scalarType, H5S_ALL, EspaceId, H5P_DEFAULT, Edata);
//     status = H5Dclose(EsetId);
   
//     EsetId = H5Dcreate(groupId, "E3", scalarType, EspaceId, H5P_DEFAULT);
//     for (int j=0; j<=J;j++)
//       for (int k=0; k<=K;k++){
//         Edata[j*(K+1) + k] = E[j][k].e3();
//       }
//     status = H5Dwrite(EsetId, scalarType, H5S_ALL, EspaceId, H5P_DEFAULT, Edata);
//     status = H5Sclose(EspaceId); 
//     status = H5Dclose(EsetId);
   
//     status = H5Gclose(groupId);
//     status = H5Fclose(fileId);

  return;
}
#endif