diagn.cpp

pic 模拟程序！面向对象

/*
====================================================================

diagn.cpp

This file contains all the arrays and functions necessary for maintaining 
and initializing the diagnostics' data.


1.1  (KC 01-6-97) Not calculation the KE of the particles here anymore,
     moved it to the push.  Add KE by species. 
2.01 (Bruhwiler 10-08-99) initialization of synchRadiationFlag added
2.02 (Bruhwiler 10-05-2000) added history plots for RMS beam quantities

====================================================================
*/

#ifdef UNIX
#include <config.h>
#endif

#include <txc_streams.h>

#include <cmath>

#include <stdlib.h>
#include <stdio.h>

#include "fields.h"
#include "grid.h"

#include "ovector.h"
#include <sptlrgn.h>
#include <ptclgrp.h>
#include "globalvars.h"
#include "oopiclist.h"
#include "diagn.h"
#include "gpdist.h"
#include "history.h"
#include "newdiag.h"


// hdf5 includes
#include "h5diagg.h"



// NPARTICLES specifies that, when the number of particles
//  of a given species exceeds this number, then XGrafix
//  plots only every 2nd particle.  It does not affect the
//  physics simulation.
#define NPARTICLES 800000
//#define HISTMAX 4096

Scalar AllDiagnostics = 1;
Scalar EnerPlots = 1;
Scalar PhaseSpacePlots = 1;

oopicList <Diagnostics> *theDiagnostics;

#ifdef MPI_VERSION
#include <mpi.h>
extern MPI_Comm XOOPIC_COMM;
extern int MPI_RANK;
#define SRB_ANNOUNCE_TAG 100
#define SRB_LINK_TAG 200
#define N_DIAGNOSTIC_TAG 300
typedef struct {
  char *name;
  int index;
  int linkedP;
} SRBdat;
#endif


//  These are the diagnostics arrays:  they're global in scope but used
//  only here and in Maintian_Diagnostics and in xgrafix.

Diagnostics::Diagnostics(SpatialRegion *thisSpace) {
  number_of_species = 1;
  Bz0=0.0;
  divderrormag = 0;
  Show_loaded_densityFlag = FALSE;
  theSpace=thisSpace;
  RegionName = theSpace->getName();
  const int jm = theSpace->getJ();
  Jm = jm;
  const int km = theSpace->getK();

  number_of_species=theSpace->get_nSpecies();
  int j,k,isp;
  B = theSpace->getBNode();
  E = theSpace->getENode();
#ifdef DEBUG_FIELDS
  intBdS = theSpace->getIntBdS();
  intEdl = theSpace->getIntEdl();
#endif
  blist = theSpace->getBoundaryList();
  BoundaryIhist = new oopicList<Ihistdiag>;
  BoundaryPFhist = new oopicList<PFhistdiag>;
  BDynamic = theSpace->getBNodeDynamic();
  I = theSpace->getI();
  rho = theSpace->getRho();

  /** Set the local pointer to oopicList<NGD>* to point
   * to the corrensponding list on the grid.
   * dad: 01/24/2001
   */
  ptrNGDList = theSpace->getPtrNGDList();

  boltzmannFlag = theSpace->getBoltzmannFlag();
  electrostaticFlag = theSpace->getElectrostaticFlag();
  synchRadiationFlag = theSpace->getSynchRadiationFlag();
  if (boltzmannFlag){
		totalRho = theSpace->getTotalRho();
	}
  loaded_density = theSpace->getloaddensity();
  Show_loaded_densityFlag = theSpace->getShowInitialDensityFlag();
  phi = theSpace->getPhi();
  divderror = theSpace->getDivDerror();
  rho_species = theSpace->get_rho_species();
  theSpecies= new SpeciesDiag[number_of_species];
  CellVolumes=theSpace->get_halfCellVolumes();

  // get memory for the grid arrays
  x1_array = new Scalar[jm+2];
  x2_array = new Scalar[km+2];

  //get memory for the diagnostics
  S_array= new Vector3 *[jm+2];
  Ue= new Scalar * [jm+2];
  Ub= new Scalar * [jm+2];
		
#ifdef DEBUG_PHI
		phi_err = new Scalar *[jm+2];
#endif

  // Set up the user-defined diagnostics
  oopicListIter<Diag> DiagIter(*thisSpace->getDiagList());
  for(DiagIter.restart();!DiagIter.Done();DiagIter++) 
	 DiagIter.current()->setDiagnostics(this);
  
  for(j=0;j<=jm;j++)
	 {
		S_array[j]= new Vector3 [km+2];
		for(k=0;k<=km;k++)
			S_array[j][k] = Vector3(0,0,0);

		Ue[j]=new Scalar[km+2];
		memset(Ue[j],0,(km+2)*sizeof(Scalar));

		Ub[j]=new Scalar[km+2];
		memset(Ub[j],0,(km+2)*sizeof(Scalar));
		
#ifdef DEBUG_PHI
		phi_err[j] = new Scalar[km+2];
#endif

	 }

  //  Get memory for the arrays that contain particle positions for xgrafix

  for(isp=0;isp<number_of_species;isp++) {
	 // search for the name corresponding to this species
	 oopicListIter<Species> siter(*theSpace->getSpeciesList());
	 for(siter.restart();!siter.Done();siter++) {
		if(isp==siter.current()->getID()) {
		  theSpecies[isp].name=siter.current()->get_name_ptr();
		  theSpecies[isp].spec = siter.current();
		  theSpecies[isp].ID = siter.current()->getID(); 
		}
	 }

	 theSpecies[isp].memory_size=NPARTICLES;
	 theSpecies[isp].nparticles=0;
	 theSpecies[isp].nparticles_plot = 0;
	 theSpecies[isp].x1positions= new Scalar[NPARTICLES];
	 theSpecies[isp].x2positions= new Scalar[NPARTICLES];
	 theSpecies[isp].x1velocities=new Scalar[NPARTICLES];
	 theSpecies[isp].x2velocities=new Scalar[NPARTICLES];
	 theSpecies[isp].x3velocities=new Scalar[NPARTICLES];
	 memset(theSpecies[isp].x1positions,0,NPARTICLES*sizeof(Scalar));
	 memset(theSpecies[isp].x2positions,0,NPARTICLES*sizeof(Scalar));
	 memset(theSpecies[isp].x1velocities,0,NPARTICLES*sizeof(Scalar));
	 memset(theSpecies[isp].x2velocities,0,NPARTICLES*sizeof(Scalar));
	 memset(theSpecies[isp].x3velocities,0,NPARTICLES*sizeof(Scalar));
  }


  number = new Scalar_History* [number_of_species];
  ke_species = new Scalar_History* [number_of_species];
  ngroups = new Scalar_History* [number_of_species];
  total_density = new Scalar_History* [number_of_species];
  Ave_KE = new Scalar_History* [number_of_species];

  HISTMAX = theSpace->get_histmax(); // retrieve time array length

  for(isp=0;isp<number_of_species;isp++) {
    number[isp]= new Scalar_History(HISTMAX,2);
    ke_species[isp] = new Scalar_History(HISTMAX,2);
    ngroups[isp] = new Scalar_History(HISTMAX,2);
    total_density[isp] = new Scalar_History(HISTMAX,2); 
    Ave_KE[isp] = new Scalar_History(HISTMAX,2);
  }
  
  /* **************************************************************
   * Need to collect scalar history of RMS beam quantities
   * for beam physics applications:  velocities, beam sizes,
   * energy spread, emittance, etc.
   * (Bruhwiler, revisions started on 10/04/2000)
   * (Bruhwiler/Dimitrov,  started on 10/09/2000)
   ****************************************************************/
  int sh;
  aveBeamSize = new Scalar_History* [2];
  rmsBeamSize = new Scalar_History* [2];
  for ( sh = 0; sh < 2; sh++ ) {
    aveBeamSize[sh]  = new Scalar_History(HISTMAX,2); 
    rmsBeamSize[sh]  = new Scalar_History(HISTMAX,2); 
  }

  // get memory for the velocity (three components in this case)
  aveVelocity = new Scalar_History* [3];
  rmsVelocity = new Scalar_History* [3];
  for ( sh = 0; sh < 3; sh++ ) {
    aveVelocity[sh] = new Scalar_History(HISTMAX,2);  
    rmsVelocity[sh] = new Scalar_History(HISTMAX,2);
  } 
  // get memory for the emittance (two components in this case)
  rmsEmittance = new Scalar_History* [2];
  for ( sh = 0; sh < 2; sh++ ) {
    rmsEmittance[sh] = new Scalar_History(HISTMAX,2);
  }  
  // get memory for the aveEnergy_eV and rmsEnergy_eV
  ///aveEnergy_eV = new Scalar_History(HISTMAX,2);
  ///rmsEnergy_eV = new Scalar_History(HISTMAX,2);

  //  Set up the grid arrays
  for(j=0;j<=jm;j++)
	 x1_array[j]=theSpace->getMKS(j,0).e1();
  for(k=0;k<=km;k++)
	 x2_array[k]=theSpace->getMKS(0,k).e2();

  //  get the size of the grid
  x1max = theSpace->getMKS(jm,km).e1();
  x2max = theSpace->getMKS(jm,km).e2();
  x1min = theSpace->getMKS(0,0).e1();
  x2min = theSpace->getMKS(0,0).e2();

  /**
   * allocate memory and initialize the 
   * x1_arrayNGD and x2_arrayNGD
   */
  jmNGD = jm;
  kmNGD = km;
  Vector2 x;  
  x1_arrayNGD = new Scalar[jm+1];
  x2_arrayNGD = new Scalar[km+1];
  for(j=0;j<jm;j++) {
    x.set_e1( static_cast<Scalar>(j) + 0.5 );
    x.set_e2( 0.5 );
    x1_arrayNGD[j] = (theSpace->getMKS(x)).e1();
  }
  for(k=0;k<km;k++) {
    x.set_e1( 0.5 );
    x.set_e2( static_cast<Scalar>(k) + 0.5 );
    x2_arrayNGD[k] = (theSpace->getMKS(x)).e2();
  }

  //set up diagnsotic history arrays for boundaries which
  //have particle diagnostics on
  oopicListIter<Boundary> nextb(*blist);
  for(nextb.restart();!nextb.Done();nextb++) { //if it's taking diagnostics of particles
     if(nextb.current()->get_particle_diag()) {
        Ihistdiag *Ihistory= new Ihistdiag(number_of_species);
        int len = nextb.current()->get_Ihist_len();
        int avg = nextb.current()->get_Ihist_avg();
        Ihistory->Ihist = new Scalar_Ave_History(len, avg);
        Ihistory->p_dist = nextb.current()->get_particle_diag();
        strcpy(Ihistory->name,nextb.current()->getBoundaryName().c_str());
        for (int i=0; i<number_of_species; i++) {
           Ihistory->Ihist_sp[i] = new Scalar_Ave_History(len,avg);
        }
        BoundaryIhist->add(Ihistory);
     }
  }
  for(nextb.restart();!nextb.Done();nextb++){ 
     if(nextb.current()->getPF()) {
        PFhistdiag *PFhistory= new PFhistdiag;
        PFhistory->PFhist = new Scalar_History(HISTMAX,2);
        PFhistory->PFlocal = new Scalar_Local_History(HISTMAX,MAX(1,HISTMAX/16));
        PFhistory->p_flux = nextb.current()->getPF();
        strcpy(PFhistory->name,nextb.current()->getBoundaryName().c_str());
        BoundaryPFhist->add(PFhistory);
     }
  }
}

Diagnostics::~Diagnostics() {

  BoundaryIhist->deleteAll();
  BoundaryPFhist->deleteAll();

  int sh;
  for ( sh = 0; sh < 2; sh++ ) {
    delete rmsEmittance[sh];
  } 
  delete [] rmsEmittance;

  for (sh = 0; sh < 3; sh++ ) {
    delete aveVelocity[sh];
    delete rmsVelocity[sh]; 
  }
  delete [] aveVelocity;
  delete [] rmsVelocity;

  for (sh = 0; sh < 2; sh++ ) {
    delete aveBeamSize[sh];
    delete rmsBeamSize[sh];
  }
  delete [] aveBeamSize;
  delete [] rmsBeamSize; 

  int isp;
  for(isp=0;isp<number_of_species;isp++) {
      delete number[isp];
      delete ke_species[isp];
      delete ngroups[isp];
      delete total_density[isp]; 
      delete Ave_KE[isp];
    }
    delete [] number; 
    delete [] ke_species;
    delete [] ngroups;
    delete [] total_density;
    delete [] Ave_KE;

  for(isp=0;isp<number_of_species;isp++) {
    delete [] theSpecies[isp].x1positions;
    delete [] theSpecies[isp].x2positions;
    delete [] theSpecies[isp].x1velocities;
    delete [] theSpecies[isp].x2velocities;
    delete [] theSpecies[isp].x3velocities;
  }
  delete [] theSpecies;

  int j;
  for(j=0;j<=Jm;j++) {
    delete [] S_array[j];
    delete [] Ue[j];
    delete [] Ub[j];		
#ifdef DEBUG_PHI
    delete [] phi_err[j];
#endif
  }
  delete [] S_array;
  delete [] Ue;
  delete [] Ub;
#ifdef DEBUG_PHI
  delete [] phi_err;
#endif  
 
  delete [] x1_array;
  delete [] x2_array;

  delete [] x1_arrayNGD;
  delete [] x2_arrayNGD;  
}

/**********************************************************************
 * a helper function to calculate averages and standard deviations
 * for an array of "random" variables distributed according to
 * the array of weight elements _Pq. It uses a two-pass algorithm.
 * DAD: Thu Oct 19 11:23:14 MDT 2000
 **********************************************************************/
void  Diagnostics::Get_Statistical_Moments(int N, Scalar* W, Scalar* Data,
					   Scalar& Ave, Scalar& Var)
{
  // first calculate the average
  Ave = 0.0;
  int i;
  for ( i = 0; i < N; i++ )
    Ave += W[i] * Data[i];
  // now calculate the variance on the second pass.
  Scalar tmp;
  Scalar sum  = 0.0;
  Scalar sum2 = 0.0;
  for ( i = 0; i < N; i++ ) {
    tmp   = Data[i] - Ave;
    sum  += W[i] * tmp;
    sum2 += W[i] * tmp * tmp;
  }

  Var = sum2 - (sum * sum);
  // check if it is positive
  if ( Var < 0.0 ) {
    if ( Var > -Scalar_EPSILON ) { // a precision problem set Var to zero 
      Var= 0.0;
    } else {
      printf(" ");
      printf("WARNING: Negative variance found in ");
	    printf("Diagnostics::Get_Statistical_Moments(...)");
	    printf("This should NEVER happen!!!");
	    printf("Setting the value of the variance to ZERO!");
    }
  }
  
  return;
}
/**********************************************************************
 * A helper function for the calculation of the averages and rms' of
 * positions, velocities, emittances, and energy for a specific 
 * species. DAD: Mon Oct 23 11:05:04 MDT 2000
 **********************************************************************/
void Diagnostics::Set_Diagnostics_Moments(int Geometry_Flag, 
			       oopicListIter<ParticleGroup>& PGIterator, 
			       Grid* RMSGrid, Scalar* _Pq, int NumWeightElements, 
			       Scalar& AveX1, Scalar& AveX2, 
			       Scalar& RMSX1, Scalar& RMSX2, 
			       Scalar* _aveU, Scalar* _rmsU, Scalar* _rmsEmit,
			       Scalar& _aveE_eV, Scalar& _rmsE_eV)
{
  int numParticles;
  int counter;
  ParticleGroup* RMSParticles = 0;  

  Vector2* position;  // pointer to particle position array