diagn.cpp

pic 模拟程序！面向对象

#if (!defined(_MSC_VER) && (  defined(__DECCXX) || defined(_CRAYT3E) || defined(__KCC) || defined(__xlC__) || defined(__ICC)))
  Scalar* X1 = new Scalar(NumWeightElements); // position arrays
  Scalar* X2 = new Scalar(NumWeightElements);
  Scalar* U1 = new Scalar(NumWeightElements); // velocity arrays
  Scalar* U2 = new Scalar(NumWeightElements);
  Scalar* U3 = new Scalar(NumWeightElements);
#else 
#if defined(__BCPLUSPLUS__) || defined(_MSC_VER) 
  Scalar* X1 = new Scalar[NumWeightElements]; // position arrays
  Scalar* X2 = new Scalar[NumWeightElements];
  Scalar* U1 = new Scalar[NumWeightElements]; // velocity arrays
  Scalar* U2 = new Scalar[NumWeightElements];
  Scalar* U3 = new Scalar[NumWeightElements];
#else 
  Scalar X1[NumWeightElements]; // position arrays
  Scalar X2[NumWeightElements];
  Scalar U1[NumWeightElements]; // velocity arrays
  Scalar U2[NumWeightElements];
  Scalar U3[NumWeightElements];
#endif  // __BCPLUSPLUS__
#endif  

  ///Scalar E_eV[NumWeightElements]; // energy array for the particles in a group in eV

  // initialize X1, X2, velocity, and energy arrays.
  Vector3* U;  // pointer to particle velocity array
  counter = 0;
  for(PGIterator.restart(); !PGIterator.Done(); PGIterator++) {
    RMSParticles = PGIterator.current();
    numParticles = RMSParticles->get_n();
    position     = RMSParticles->get_x();
    U            = RMSParticles->get_u();
    for (int j=0; j<numParticles; j++) {
      X1[counter]   = RMSGrid->getMKS( position[j] ).e1();
      X2[counter]   = RMSGrid->getMKS( position[j] ).e2();
      U1[counter]   = U[j].e1();
      U2[counter]   = U[j].e2();
      U3[counter]   = U[j].e3();
      ///E_eV[counter] = RMSParticles->energy_eV(j);
      counter++;
    }
  }
  
  // get the mean and standard deviation values for the positions
  Scalar VarX1, VarX2;
  Get_Statistical_Moments(NumWeightElements, _Pq, X1, AveX1, VarX1);
  Get_Statistical_Moments(NumWeightElements, _Pq, X2, AveX2, VarX2);
  RMSX1 = sqrt( VarX1 );
  if ( Geometry_Flag == 0 ) { // cylindrical geometry according to "grid.h"   
    // average of X2 == r must be zero since it is the average
    // of the vector r. (Bruhwiler/Dimitrov, 10/19/2000)
    VarX2 += AveX2*AveX2; 
    RMSX2 = sqrt( VarX2 );
    AveX2 = 0.; 
  } else if ( Geometry_Flag == 1 ) { // xy geometry
    RMSX2 = sqrt( VarX2 );
  }
  
  // get the mean and standard deviation values for the velocities
  Scalar VarU[3];
  Get_Statistical_Moments(NumWeightElements, _Pq, U1, _aveU[0], VarU[0]);
  Get_Statistical_Moments(NumWeightElements, _Pq, U2, _aveU[1], VarU[1]);
  Get_Statistical_Moments(NumWeightElements, _Pq, U3, _aveU[2], VarU[2]);
  for ( int j = 0; j < 3; j++ ) 
    _rmsU[j] = sqrt(VarU[j]);

  // calculate the rms Emittance for X1 and X2;  
  if ( Geometry_Flag == 0 ) { // cylindrical geometry according to "grid.h"   
    // by default for cylindrical geometry the beam axis is X1
    // this will be calculated only if fabs(_aveU[1]/_aveU[0]) < 1.0e-3
    // and fabs(_aveU[2]/_aveU[0]) < 1.0e-3 which will be use as a condition 
    // for a beam, otherwise the emittance will be set to the negative
    // value of 0.0
    if ( (fabs(_aveU[1]/_aveU[0]) < 1.0e-3) && (fabs(_aveU[2]/_aveU[0]) < 1.0e-3) ) {
      Scalar ave_XU[2];
      int i;
      for ( i = 0; i < 2; i++ ) {
	ave_XU[i] = 0.0;
      }
      for ( i = 0; i < NumWeightElements; i++ ) {
	ave_XU[0] += _Pq[i] * (X1[i] - AveX1) * (U1[i] - _aveU[0]);
	ave_XU[1] += _Pq[i] * (X2[i] - AveX2) * (U2[i] - _aveU[1]);
      }
      
      _rmsEmit[0] = VarX1 * VarU[0] - ave_XU[0] * ave_XU[0];
      _rmsEmit[1] = VarX2 * VarU[1] - ave_XU[1] * ave_XU[1];
      
      for ( i = 0; i < 2; i++ ) {
	if ( _rmsEmit[i] < 0.0 ) {
	  if ( _rmsEmit[i] > -Scalar_EPSILON ) { /* a precision problem,  
						    set _rmsEmit[i] to zero */
	    _rmsEmit[i] = 0.0;
	  } else {
          TXSTRSTD::cout << "WARNING: Negative value for the square of the " 
	       << "Emittance E[" << i+1 << "] found in " << endl 
	       << "Diagnostics::Set_RMS_Cylindrical(...)" << endl
	       << "This should NEVER happen!!! Now setting the value of" << endl
	       << "the emittance to ZERO!" << endl;
	  }
	} else {
	  _rmsEmit[i] = sqrt(_rmsEmit[i]) / fabs( _aveU[0]);
	}
      }
    } else {
      for ( int i = 0; i < 2; i++ ) 
	_rmsEmit[i] = 0.0;
    }
  } else if ( Geometry_Flag == 1 ) { // xy geometry
    // for cartesian geometry we have to determine the 
    // direction of the beam, if we have a beam at all. 
    static int BeamDirectionIndex = -1; // assume initially there is no beam
    int AveUmaxIndex;
    int ib = 0;
    ib = (fabs(_aveU[0]) > fabs(_aveU[1]) ) ? 0 : 1;
    AveUmaxIndex = (fabs(_aveU[ib]) > fabs(_aveU[2])) ? ib : 2;
    // is it a beam?
    if ( AveUmaxIndex != ib ) {
      if ( fabs(_aveU[ib]/_aveU[AveUmaxIndex]) < 1.0e-3 )
	BeamDirectionIndex = AveUmaxIndex;
    } else {
      if ( AveUmaxIndex == 0 )
	ib = (fabs(_aveU[1]) > fabs(_aveU[2]) ) ? 1 : 2;
      else if ( AveUmaxIndex == 1 ) 
	ib = (fabs(_aveU[0]) > fabs(_aveU[2]) ) ? 0 : 2;
      if ( fabs(_aveU[ib]/_aveU[AveUmaxIndex]) < 1.0e-3 )
	BeamDirectionIndex = AveUmaxIndex;
    } 
    
    int i;
    if ( BeamDirectionIndex == -1 ) { // no beam
      for ( i = 0; i < 2; i++ ) 
	_rmsEmit[i] = 0.0; // set it ZERO for no beam
    } else { 
      // there is a beam
      // calculate the rms Emittance for X1 and X2;
      Scalar ave_XU[2];
      for ( i = 0; i < 2; i++ ) {
	ave_XU[i] = 0.0;
      }
      for ( i = 0; i < NumWeightElements; i++ ) {
	ave_XU[0] += _Pq[i] * (X1[i] - AveX1) * (U1[i] - _aveU[0]);
	ave_XU[1] += _Pq[i] * (X2[i] - AveX2) * (U2[i] - _aveU[1]);
      }
      
      _rmsEmit[0] = VarX1 * VarU[0] - ave_XU[0] * ave_XU[0];
      _rmsEmit[1] = VarX2 * VarU[1] - ave_XU[1] * ave_XU[1];
      for ( i = 0; i < 2; i++ ) {
	if ( _rmsEmit[i] < 0.0 ) {
	  if ( _rmsEmit[i] > -Scalar_EPSILON ) { /* a precision problem,  
						    set _rmsEmit[i] to zero */
	    _rmsEmit[i] = 0.0;
	  } else {
	  TXSTRSTD::cout << "WARNING: Negative value for the square of the " 
	       << "Emittance E[" << i+1 << "] found in " << endl 
	       << "Diagnostics::Set_RMS_Cylindrical(...)" << endl
	       << "This should NEVER happen!!! Now setting the value of" << endl
	       << "the emittance to ZERO!" << endl;
	  }
	} else {
	  _rmsEmit[i] = sqrt(_rmsEmit[i]) / fabs( _aveU[0]);
	}
      }
    }
  }

  // get the mean and standard deviation values for the Energy
  ///Scalar varE;
  ///Get_Statistical_Moments(NumWeightElements, _Pq, E_eV, _aveE_eV, varE);
  ///_rmsE_eV = sqrt(varE);

  // Clean up the locally allocated memory:
#if defined(__DECCXX) || defined(__BCPLUSPLUS__) || defined(_CRAYT3E) || defined(__KCC) || defined(_MSC_VER) || defined(__xlC__)
  delete [] X1;
  delete [] X2;
  delete [] U1;
  delete [] U2;
  delete [] U3;
#endif
  
  return;
} 

/************************************************************************/
/* time history accumulator; calculates and stores all history values,  */
/* and performs combing on history values when low on memory            */
/************************************************************************/
void Diagnostics::history()
{
	register int isp, k, j;

	/*****************************************/
	/*  Fixing the diagnostic arrays         */

	simulation_time=theSpace->getTime();
	int jm=theSpace->getJ();
	int km=theSpace->getK();
	oopicListIter<Ihistdiag> nextdiag(*BoundaryIhist);
	oopicListIter<PFhistdiag> nextdiag2(*BoundaryPFhist);

	DiagList *dlist = theSpace->getDiagList();
	oopicListIter<Diag> nextd(*dlist);

	// ------Updating newDiag-----//
	for(nextd.restart();!nextd.Done(); nextd++) 
		nextd.current()->MaintainDiag((Scalar)simulation_time);
	//--------------------------//

	// update the history objects

	Uktot = 0;
	for(isp=0;isp<number_of_species;isp++) {
		number[isp]->add(theSpecies[isp].nparticles,simulation_time);
		ngroups[isp]->add(theSpecies[isp].ngroups,simulation_time);
		ke_species[isp]->add(theSpecies[isp].spec->getKineticEnergy(),simulation_time);
		Uktot += theSpecies[isp].spec->getKineticEnergy();
	}
	energy_e.add(Uetot,simulation_time);
	energy_b.add(Ubtot,simulation_time);
	energy_k.add(Uktot,simulation_time);
	energy_all.add(Uetot+Ubtot+Uktot,simulation_time);
	divderrorhis.add(divderrormag,simulation_time);

	for(nextdiag2.restart();!nextdiag2.Done();nextdiag2++)
		{ //if it's taking diagnostics of Poyting Flux
			Scalar p_flux = *nextdiag2.current()->p_flux;
			nextdiag2.current()->PFhist->add(p_flux,simulation_time);
			nextdiag2.current()->PFlocal->add(p_flux,simulation_time);
		}
	for(nextdiag.restart();!nextdiag.Done();nextdiag++) { 
	  //if it's taking diagnostics of particles
	  Scalar Ir = nextdiag.current()->p_dist->getdQ()/dt;
	  nextdiag.current()->Ihist->add(Ir,simulation_time);
	  for (isp=0; isp<number_of_species; isp++) {
	    Ir = nextdiag.current()->p_dist->getdQ(isp)/dt;
	    nextdiag.current()->Ihist_sp[isp]->add(Ir,simulation_time);
	  }
	}

	// Calculation of the total density
	Scalar total_dens = 0.0;
	Scalar icharge = 0.0;
	for(isp=0;isp<number_of_species;isp++) {
		total_dens = 0.0;
		icharge = 1/theSpecies[isp].spec->get_q(); 
		for (j=0; j<=jm; j++)
			for (k=0; k<=km; k++)
				total_dens += CellVolumes[j][k]*rho_species[isp][j][k];

		total_dens *= icharge;

		total_density[isp]->add(total_dens,simulation_time);
		if (total_dens)
			Ave_KE[isp]->add(theSpecies[isp].spec->getKineticEnergy()/(fabs(ELECTRON_CHARGE)*total_dens), simulation_time);
		else
			Ave_KE[isp]->add(0.,simulation_time);
	}

	/***************************begin********************************
	 * Calculation of the rms beam size (presumably for a beam)
	 * (Bruhwiler/Dimitrov, 10/09/2000)
	 ****************************************************************/
	// This is done for species number iSpecesRMS, assuming that the
	// rmsBeamSizeFlag has been set to 1 (default value 0) for at
	// least one species.  If the flag has been set for two different
	// species, only the first one will get plotted.
	// We set the initial value to -1, so that nothing is plotted,
	// if the rmsBeamSizeFlag has not been set to 1 for any species.
	int iSpeciesRMS = -1;
	
	// We will loop over all species in the simulation.
	// For each species, we will check whether the
	// rmsDiagnosticsFlag has been set to 1 (default value is 0).
	oopicList<ParticleGroup> *pgList;
	oopicListIter<ParticleGroup> pgIterator;
	ParticleGroup* tmpGroup;
	Species* tmpSpecies;
	int tmpFlag;
	BOOL Species_Flag = false;
	isp = 0; 
	while ( (isp < number_of_species ) && (!Species_Flag) ) {
	  oopicListIter<ParticleGroup> 
	    pgscan(*theSpace->getParticleGroupListArray()[isp]);
	  pgscan.restart();
	  if ( pgscan.isEmpty() ) {
	    tmpGroup = pgscan.current();
	    tmpSpecies = tmpGroup -> get_species();
	    tmpFlag = tmpSpecies -> get_rmsDiagnosticsFlag();
	    
	    if ( tmpFlag ) {
	      iSpeciesRMS = isp;
	      pgList = theSpace->getParticleGroupListArray()[iSpeciesRMS];
	      pgIterator = oopicListIter<ParticleGroup>(*pgList);
	      Species_Flag = true;
	    }
	  }
	  isp++;
	}

	if (iSpeciesRMS != -1) {
	  // initialize a Scalar array with the 
	  // normalized the weight elements (distribution density)
	  // 
	  // (i) count the number of weight elements
	  ParticleGroup* RMSParticles = 0;
	  int numParticles; 
	  int numWeightElements = 0;
	  for(pgIterator.restart(); !pgIterator.Done(); pgIterator++) {
	    RMSParticles = pgIterator.current();
	    numWeightElements += RMSParticles->get_n();
	  }
	  // (ii) allocate memory for the weight elements
	  Scalar* Pq = new Scalar[numWeightElements];
	  // (iii) initialize Pq and the normalization constant
	  numWeightElements = 0;
	  Scalar sum_q   = 0.;
	  for(pgIterator.restart(); !pgIterator.Done(); pgIterator++) {
	    RMSParticles = pgIterator.current();
	    numParticles = RMSParticles->get_n();
	    for (int j=0; j<numParticles; j++) {
	      Pq[numWeightElements] = fabs( RMSParticles->get_q(j) );
	      sum_q   += Pq[numWeightElements];
	      numWeightElements++;
	    }
	  }
	  // (iv) normalize Pq
	  for (int j=0; j< numWeightElements; j++) {
	      Pq[j]  /= sum_q;
	  }

	  // the Grid knows what geometry it is defined for via the variable
	  // "int geometryFlag" which is a data member of the Grid class. 
	  // By default this variable is set to 0 == cylindrical geometry.
	  // It can be set to cartesian geometry via adding the line
	  // "Geometry = 1" in the input file when the Grid parameters are 
	  // specified. DAD::Wed Oct 11 14:59:01 MDT 2000
	  Grid* rmsGrid = theSpace->get_grid();
	  int tmp_Geometry_Flag = rmsGrid->query_geometry();
	  
	  // Next, we use helper functions to calculate the Ave & RMS values.
	  Scalar aveX1;
	  Scalar aveX2;
	  Scalar rmsX1;
	  Scalar rmsX2;
	  Scalar aveU[3];
	  Scalar rmsU[3];
	  Scalar rmsEmit[2]; // for the Emittance
	  Scalar aveE_eV;    // for the average and rms Energy in eV
	  Scalar rmsE_eV;
	  
	  if ( (tmp_Geometry_Flag != 0) && (tmp_Geometry_Flag != 1) ) { 
	    printf("WARNING: Undefined geometry flag detected in diagn.cpp");
	  } else {
	    Set_Diagnostics_Moments(tmp_Geometry_Flag, pgIterator, rmsGrid, Pq, 
				    numWeightElements, aveX1, aveX2, rmsX1, rmsX2, 
				    aveU, rmsU, rmsEmit, aveE_eV, rmsE_eV);   
	  }
          
	  aveBeamSize[0]->add(aveX1,simulation_time);
	  aveBeamSize[1]->add(aveX2,simulation_time);	  
	  rmsBeamSize[0]->add(rmsX1,simulation_time);
	  rmsBeamSize[1]->add(rmsX2,simulation_time);

	  aveVelocity[0]->add(aveU[0],simulation_time);
	  aveVelocity[1]->add(aveU[1],simulation_time);
	  aveVelocity[2]->add(aveU[2],simulation_time);
	  
	  rmsVelocity[0]->add(rmsU[0],simulation_time);
	  rmsVelocity[1]->add(rmsU[1],simulation_time);
	  rmsVelocity[2]->add(rmsU[2],simulation_time);

	  rmsEmittance[0]->add(rmsEmit[0],simulation_time);
	  rmsEmittance[1]->add(rmsEmit[1],simulation_time);

	  ///aveEnergy_eV->add(aveE_eV,simulation_time);
	  ///rmsEnergy_eV->add(rmsE_eV,simulation_time);
	  
	  delete[] Pq;
	}
	/****************************************************************
	 * (end of changes by Bruhwiler/Dimitrov, 10/09/2000)
	 ****************************************************************/

	return;
}

/* another function for testing. */
void Diagnostics::compute_div_derror_magnitude(Scalar **divderror,int jm,int km) {
	int j,k;
	Scalar rhomagave;
	divderrormag=0;
	rhomagave = 0;
	for(j=1;j<jm-1;j++)
		for(k=1;k<km-1;k++) {
			rhomagave+=fabs(rho[j][k]);
			divderrormag+=fabs(divderror[j][k]);
		}
	if(rhomagave!=0) divderrormag/=rhomagave;
}


/*  UpdatePreDiagnostics does processing on raw physics data such
	which should be done before a physics iteration so that the user
	gets fields + particle positions which are reasonably synchronized.
 */

void Diagnostics::UpdatePreDiagnostics()
{  
	int i,j;

#ifndef PHI_TEST

	if(AllDiagnostics) {  //collect diagnostics at all?
		 
		if(PhaseSpacePlots!=0.0) {
			int isp;
			int skip;
			//  maintain the phase space diagnostics.  
			for(isp = 0; isp < number_of_species; isp++) {
				SpeciesDiag *thisSpc= &theSpecies[isp];
				thisSpc->nparticles=0;
				int pgnumparticles;
				oopicListIter<ParticleGroup> 
					pgscan(*theSpace->getParticleGroupListArray()[isp]);
				i=0;  //place in phase space diagnostic


				// First count the particles, store in thisSpc=theSpecies[isp]