ptclgrp.cpp

pic 模拟程序！面向对象

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

File:		ptclgrp.cpp

Purpose:	Implementation of a ParticleGroup.

Version:	$Id: ptclgrp.cpp,v 1.32 2005/07/27 11:47:16 rmtrines Exp $

Revision History
0.99	(JohnV 01-07-94) Separate from pic.h.
0.991	(JohnV 03-23-94) Streamline includes, remove pic.h.
0.992	(JohnV, KC 05-02-94) Fix Brilloiun beam problem (was neglecting
      centrifugal force).
0.993	(JohnV 09-13-94) Add variable np2c support.
0.994 (PeterM 09-04-95) Added a method to reduce #particles, increase
      the weight of particles.
2.001 (JohnV 11-10) Improve increaseWeight(), handle odd n. 
2.02  (JohnV 12-12-97) Modified increaseParticleWeight() to retain temporal
      order of particles within groups -> handles default.inp properly.
2.03  (KC 1-9)  Added Kinetic Energy diagnostic to advance.
2.04  (Bruhwiler 9/29/99) set bool synchRadiationOn=FALSE in constructor
2.05  (Bruhwiler 11/08/99) now moving window only checks for particles
      leaving the left boundary of the simulation
2.0?  (Cary 23 Jan 99) Noted in a comment that this coding is only for
	right moving windows.
CVS-1.16.2.20 (JohnV 15 Mar 00) Fixed bug in np2cFactor (getting ignored)
2.1   (JohnV 10Jul2003) fixed recomputation of gamma before x update

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

#include "ptclgrp.h"
#include "particle.h"
#include "fields.h"
#include "dump.h"

#include <txc_streams.h>
using namespace std;

ParticleGroup::ParticleGroup(int _maxN, Species* s, Scalar _np2c0, 
	BOOL vary_np2c)
{
	maxN = _maxN;
	if (vary_np2c)
		qArray = new Scalar[maxN];
	else qArray = 0;	
	np2c0 = _np2c0;
	species = s;
	q_over_m = species->get_q_over_m();
	q = np2c0*species->get_q();
//	speciesSub = species->get_speciesSub();
	x = new Vector2[maxN];
	u = new Vector3[maxN];
	n = 0;

   synchRadiationOn = FALSE;  // by default, synchrotron radiation is off
}

ParticleGroup::~ParticleGroup()
{
	n = maxN = 0;								// deallocate memory
	delete[] x;
	delete[] u;
	if (qArray)
	  delete[] qArray;
}

//--------------------------------------------------------------------
//	This add uses nominal np2c0, q, q_over_m

BOOL ParticleGroup::addMKS(Vector2 _x, Vector3 _v)
{
	if (n < maxN)
	{
		x[n] = _x;
		Vector3	beta = iSPEED_OF_LIGHT*_v;
		u[n] = _v/sqrt(1 - beta*beta);	//	u = gamma*v
		if (qArray) qArray[n] = q;
		n++;
		return TRUE;
	}
	return FALSE;
}

//--------------------------------------------------------------------
//	Add a group of particles.  Not well tested yet.

BOOL ParticleGroup::add(ParticleGroup& p)
{
	if (p.q_over_m == q_over_m && p.np2c0 == np2c0)
		while (--p.n >= 0)
		{
		  if (p.qArray){
			 if (!add(x[p.n], u[p.n], p.get_np2c(p.n)))
				break;
		  }
		  else if (p.get_q(p.n) == q){
			 if (!add(x[p.n], u[p.n]))	//	did add fail?
				break;						//	bail out
		  }
		}
	p.n++;										//	restore proper value
	if (p.n) return FALSE;
	else return TRUE;     					//	caller responsible for destructor?
}

//-------------------------------------------------------------------
// duplicates particles in group N times, dividing the particle weight appropriately

/*
ParticleGroup**     ParticleGroup::duplicateParticles(int N)
{
  int i, j, ntemp;
  BOOL varweight;
  Vector2 xN;
  Vector3 uN;
  Scalar np2cfactor;
  Scalar Nfactor;
  ParticleGroup **pg;
  ParticleGroup *p;

  if(N < 1) { return NULL; }
  Nfactor = 1.0/((Scalar)(N+1));
  ntemp = n;
  varweight = isVariableWeight();

  //  fprintf(stderr, 
  //	  "\nduplicateParticles %d times, Nfactor=%g\n", N, Nfactor);
  //  cerr << "\nVariable weight: " << varweight << ", np2c0 =" << np2c0;

  pg = new ParticleGroup*[N];

  if(ntemp > 0)
    {
      if(varweight)
	{
	  np2cfactor = np2c0;
	}
      else
	{
	  np2cfactor = np2c0*Nfactor;
	}
    }
  else
    {
      return NULL;
    }

  //  fprintf(stderr, "\nnp2cfactor=%g", np2cfactor);

  for(i=0; i < N; i++)
    {
      pg[i] = new ParticleGroup(maxN, species, np2cfactor, isVariableWeight());
    }
  
  for(j=0; j < ntemp; j++)
    {
      xN = x[j];
      uN = u[j];
      np2cfactor = get_np2c(j)*Nfactor;
      if(varweight) 
	{ 
	  qArray[j] = q*np2cfactor/np2c0;
	}
      
      for(i=0; i < N; i++)
	{
	  p = pg[i];

	  if(varweight)
	    {
	      if(!p->add(xN, uN, np2cfactor)) 
		{ 
		  fprintf(stderr, "\nduplicateParticles: error adding particles");
		  return NULL; 
		}
	    }
	  else
	    {
	      if(!p->add(xN, uN, -1.)) 
		{ 
		  fprintf(stderr, "\nduplicateParticles: error adding particles");
		  return NULL; 
		}
	    }	  
	}
    }
  if(!varweight)
    {
      np2c0 *= Nfactor;
      q = np2c0*species->get_q();
      //      fprintf(stderr, "\nSetting base particle weight to %g\n", np2c0);
    }
  return pg;
}
*/

ParticleGroup**     ParticleGroup::duplicateParticles(int N)
{
  int i, j, ntemp;
  BOOL varweight;
  Vector2 xN;
  Vector3 uN;
  Scalar np2cfactor;
  Scalar Nfactor;
  ParticleGroup **pg;
  ParticleGroup *p;

  if(N < 1) { return NULL; }
  Nfactor = 1.0/((Scalar)(N+1));
  ntemp = n;
  varweight = isVariableWeight();

  //  fprintf(stderr, 
  //	  "\nduplicateParticles %d times, Nfactor=%g\n", N, Nfactor);
  //  cerr << "\nVariable weight: " << varweight << ", np2c0 =" << np2c0;

  pg = new ParticleGroup*[N];

  if(ntemp > 0)
    {
      np2cfactor = np2c0;
    }
  else
    {
      return NULL;
    }

  //  fprintf(stderr, "\nnp2cfactor=%g", np2cfactor);

  for(i=0; i < N; i++)
    {
      pg[i] = new ParticleGroup(maxN, species, np2cfactor, isVariableWeight());
    }
  
  for(j=0; j < ntemp; j++)
    {
      xN = x[j];
      uN = u[j];
      np2cfactor = get_np2c(j);
      if(varweight) 
	{ 
	  qArray[j] = q*np2cfactor/np2c0;
	}
      else
	{
	  np2cfactor = -1.;
	}
      
      for(i=0; i < N; i++)
	{
	  p = pg[i];

	  if(!p->add(xN, uN, -1.)) 
	    { 
	      fprintf(stderr, "\nduplicateParticles: error adding particles");
	      return NULL; 
	    }
	}
    }

  return pg;
}

//--------------------------------------------------------------------
//	Add a new particle to this group, possibly of variable weight.
//	Note that uNew is gamma*v, and xNew is in grid coords.

BOOL ParticleGroup::add(const Vector2& xNew, const Vector3& uNew, Scalar np2cNew)
{
	if (n < maxN)
	{
		if (np2cNew > 0)
			if (qArray)	qArray[n] = q*np2cNew/np2c0;
			else return FALSE;
		x[n] = xNew;
		u[n] = uNew;
		n++;
		return TRUE;
	}
	return FALSE;
}

//--------------------------------------------------------------------
//	Add a Particle to this group.

BOOL ParticleGroup::add(Particle* p)
{
//	if (p->get_q_over_m() == q_over_m)
	if (p->get_speciesID() == get_speciesID() && n < maxN)
		if (qArray && p->isVariableWeight())
			//		if (qArray)
			return add(p->get_x(), p->get_u(), p->get_np2c());
		else if (p->get_np2c() == np2c0)
			return add(p->get_x(), p->get_u());
	return FALSE;
}

//
// we need to be able to remove particles from the particle group (in the
// tunneling ionization of multi electron atoms). Note the implementation
// of this function: a particle with index particleIndex is deleted and
// the last particle in the group is placed at its index.    
// Returns TRUE if the particle group becomes empty (at deletion of 
// the last particle in the group), otherwise it resturns FALSE. dad 06/29/01.
//
BOOL ParticleGroup::remove(int particleIndex) {
  //
  // is particleIndex in the allowable range
  // 
  if ( particleIndex < 0 || particleIndex > (n-1) ) {
    cerr << endl << endl << endl
         << "An out of range particleIndex = " << particleIndex << endl
         << "passed to BOOL ParticleGroup::remove(int particleIndex)." << endl
         << "This index must be in [0, " << n-1 << "]" << endl
         << "Exiting to the system..." << endl << endl << endl;
  }
  
  if ( n == 1 ) {
    // 
    // we are deleting the only macroparticle in the group
    // => return false to prompt the deletion of the group
    // 
    n--;
    return TRUE;
  } else if ( particleIndex == (n-1) ) {
    //
    // deleting the last macroparticle in a group with more
    // than one macroparticles => need only to decrement the
    // number of macroparticles in the group
    //  
    n--;
    return FALSE;
  } else {
    // 
    // deleting any other particle but the last one
    // 
    x[particleIndex] = x[n-1];
    u[particleIndex] = u[n-1];
    if ( qArray ) 
      qArray[particleIndex] = qArray[n-1];
    n--;
    return FALSE;
  } 
}

//--------------------------------------------------------------------
//	Advance particle equations of motion one timestep.  This includes
//	weighting fields to the particle, solving the equation of motion,
//	and accumulating the particle current on the grid.  Boris mover
//	with u = gamma*v, v in MKS units; see B&L for variables.

void ParticleGroup::advance(Fields& fields, Scalar dt)
{
	Scalar	f = 0.5*dt*q_over_m;
	Scalar	q_dt = q/dt;
	Vector3	uPrime;
	Vector3	a;
	Vector3	t;
	Vector3	s;
	Vector2*	x = ParticleGroup::x;
	Vector3*	u = ParticleGroup::u;
	Vector3	dxMKS;
	Boundary*	bPtr;
	Grid&	grid = *fields.get_grid();
	Scalar igamma;
	Scalar localEnergy = 0;

	for (int i=0; i<n; i++, x++, u++)
	{
		if (qArray) q_dt = qArray[i]/dt;	//	variable weighted particles
		a = f*fields.E(*x);
//
		*u += a; // half acceleration
		igamma = 1/ParticleGroup::gamma(*u);
		t = (f*igamma)*fields.B(*x);
		uPrime = *u + u->cross(t);
		s = 2*t/(1 + t*t);
		*u += uPrime.cross(s); // rotation
		*u += a; // half acceleration
		igamma = 1/ParticleGroup::gamma(*u); // must recompute gamma
		dxMKS = grid.differentialMove(fields.getMKS(*x), *u*igamma, *u, dt);
		if ((bPtr = fields.translateAccumulate(*x, dxMKS, q_dt)))
		{
			//	send this particle to boundary
			bPtr->collect(*(new Particle(*x, *u, species, get_np2c(i),
				(BOOL)(qArray!=0))), dxMKS);
			//	Move last particle into this slot and advance it.
			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];
			i--;									//	Set index to do i again.
			x--;									//	Setup to redo the array element
			u--;                          // it gets inc in for loop
		}
		else{
			if (qArray)
				localEnergy += qArray[i]*(1/igamma-1);
			else
				localEnergy += 1/igamma-1;
		}
#ifdef DEBUG
		if(strncmp(HOSTTYPE, "alpha", 5) != 0)
		  {
		    if(((*x).e1()==grid.getJ() || (*x).e1()==0 || (*x).e2()==grid.getK()
			|| (*x).e2()==0)&&!bPtr) {
		      printf("Particle escaped.\n");
		      int *crash = 0;	
		      *crash = 10;
		    }
		  }
#endif //DEBUG
	}
	if (qArray)
		species->addKineticEnergy(SPEED_OF_LIGHT_SQ*localEnergy/q_over_m);
	else
		species->addKineticEnergy(SPEED_OF_LIGHT_SQ*localEnergy*get_m());
}

// this function combs the particles, doubling weights
// and killing every second particle. Note that the particles must
// remain time-ordered so that multiple packings do not leave gaps.

void ParticleGroup::increaseParticleWeight()
{
	int new_n = (int) (n*0.5 + frand()); // keep odd particles 50% time
	Scalar ratio = 2; // = ((Scalar) n)/new_n;
	if (qArray) qArray[0] *= ratio;
	for (int i = 1; i < new_n; i++)
	{
		x[i]=x[2*i];
		u[i]=u[2*i];
		if(qArray) qArray[i] = qArray[2*i]*ratio;
	}
	n = new_n;	// halve the num of particles
	np2c0 *= ratio;	// double the weights
	q *= ratio;		// double the charge
}

void ParticleGroup::fill(int index)
{
	n--;
	if (index == n) return;				//	last particle in array?
	x[index] = x[n];		//	move last particle into open slot
	u[index] = u[n];
	if (qArray) qArray[index] = qArray[n];
}
#ifdef HAVE_HDF5
int ParticleGroup::dumpH5(dumpHDF5 &dumpObj, Grid *grid,string groupName) 
{

  //cerr << "entering ptclgrp::dumpH5(dumpObj)\n";