sptlrgn.cpp

pic 模拟程序！面向对象

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

SPTLRGN.CPP

Purpose:	functions of the SpatialRegion class.  The SpatialRegion
		is the fundamental building block of the PlasmaDevice to be
		simulated.

Revision/Author/Date
0.1	(JohnV, 03-15-93)	Proto code.
0.9	(JohnV, 04-30-93)	Initial draft with particles.
0.91	(NTG, JohnV, 05-05-93) Add EmittingBoundary, ActiveBoundary, class
		library use for list management and iterators, KinematicalObject.
0.93	(JohnV, NTG 05-12-93) Fix Windows large model Gen. Prot. Fault due
		to reading outside arrays.
0.95	(JohnV, 06-29-93) Add Boundary::setPassives() to set up passive
		boundary conditions (Fields::iC and ::iL), add Fields::set_iCx()
		to support above.
0.96	(JohnV 08-23-93) Add new Boundary subtypes: BeamEmitter,
		CylindricalAxis, FieldEmitter, Port.  Modifications to
		Maxwellian: stores all parameters in MKS.
0.961	(PeterM, JohnV, 08-24-93) Replace Borland list containers with
		vanilla C++ template lists in oopiclist.h
0.963	(JohnV, 09-24-93)	Add Particle object for moving particle quantities
		around (emission, collection, etc.)
0.964	(JohnV, 11-22-93) Make agreed-upon changes in names: Vector->Vector3,
		Boundaries->BoundaryList, Particles->ParticleList, ParticleGroups
		->ParticleGroupList
0.965	(JohnV, 02-01-94) Move magnetic field advance to Fields::advanceB(),
		move boundary code, particle code to respective files.
0.966	(JohnV, 02-07-94)	Fixed striation in vr vs. z phase space seen
		for beam spreading in a simple cylinder.  Problem was weighting
		error in Grid::interpolateBilinear().
0.967	(JohnV 03-03-94) added Fields::setBz0() and ::epsilonR.
0.968	(JohnV 03-11-94) Fix divide by zero for non-translating particles
		with finite spin in Fields::translateAccumulate().  Also ensure all
		particles that touch boundaries get collected in Grid::translate()
		by changing </> to <=/>= for check conditions.
0.969	(JohnV 03-23-94) Restructure Grid, Fields, and SpatialRegion to
		an association rather than inheritance.
0.970	(GMU 11-16-94) Added TaskYield() call.
0.971 (JohnV 01-26-95) Add code to handle ParticleGroup array and Species.
0.98  (KeithC 08-25-95) Lack of precision caused Bz0 to 'leak` into Etheta 
      and Br.  Added BNodeStatic to separate the Static part from the B
		fields that that Fields  updates.  Changed Fields::setBz0() to 
		Fields::setBNodeStatic.
0.981 (PeterM 09-04-95) Added a method to reduce #particles, increase
      the weight of particles.
1.0	(JohnV 01-18-96) Added subcycling stuff from CooperD; rewrote advance().
1.1   (JohnV 01-21-97) Added supercycling.
1.11  (JohnV 10-01-97) Improved supercycling to work with mcc package.
2.01  (JohnV 11-10-97) CountParticles() deletes empty ParticleGroup's.
2.02  (JohnV 12-12-97) Modified increaseParticleWeight().
2.03  (JohnV 03-31-98) Implement caching scheme for nParticles
2.04  (Bruhwiler 10-08-99) Added getSynchRadiationFlag() check to addParticleList()
2.0?  (Cary 22 Jan 00) Added methods for setting the shift delay time and
	streamlined the shifting code.
CVS-1.53.2.7 (Cary, 22 Jan 00) reorganized calls for shifting in
	anticipation of making MPI work.
CVS-1.53.2.8 (Cary, 22 Jan 00) Added code to determine which boundary is
	involved.  The if is a SPATIAL_REGION_BOUNDARY, call passParticles().
	Similar coding needed for receiving the shifted particles and fields.
	I am not sure that this works except for right and down shifting
	windows.
CVS-1.53.2.12 (Cary 27 Jan 00) ShiftRegion changed to shiftRegion.
CVS-1.53.2.13 (kgl  28 Jan 00) xgmini.h is not needed, only theDumpFile decl.)
CVS-1.53.2.14 (JohnV 15 Mar 00) Fixed bug in np2cFactor (getting ignored)
CVS-1.57      (kgl  4  Aug 00) dropfile size (50) is too small, increased
                               to 80.  Caused coredump for long dir.
                               Changed to 180 by kgl and dad to accomodate
                               long file names. 04/08/02.

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

//I/O stuff

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

#include <txc_streams.h>
using namespace std;

// Standard includes
#include <cassert>

// MPI stuff
#ifdef MPI_VERSION
#include <unistd.h>
#include <mpi.h>

extern int MPI_RANK;

#ifdef MPI_DEBUG
extern std::ofstream dfout;
#endif  // MPI_DEBUG

#endif  // MPI_VERSION


#include "sptlrgn.h"
#include "ptclgrp.h"
#include "ptclgpes.h"
#include "ptclgpsr.h"
#include "ptclgpnr.h"
#include "ptclgpib.h"
#include "mcc.h"
#include "periodic.h"
#include "load.h"
// #include "xgmini.h"
#ifdef OOPICPRO

#if defined(_MSC_VER)
extern "C" char theDumpFile[180];
#else
extern char theDumpFile[180];
#endif

#else

#if defined(_MSC_VER)
extern "C" char* const volatile theDumpFile;
#else
extern char* const volatile theDumpFile;
#endif

#endif
// for dumping diagnostics

#include "diagn.h"
extern oopicList<Diagnostics> *theDiagnostics;



//--------------------------------------------------------------------
//	Build SpatialRegion, along with inherited Fields object

extern void taskYield();    //Yield execution to concurrent processes

SpatialRegion::SpatialRegion(Fields* f, BoundaryList* b, 
	ParticleGroupList** p, SpeciesList* s, MCC* _mcc, 
                             MCTI* _mcti, Scalar _dt) {

  fields = f;
  grid = fields->get_grid();
  boundaryList = b;
  particleGroupList = p;

 

  speciesList = s;
  nSpecies = speciesList->nItems();
  speciesArray = new Species*[nSpecies];
  nParticles = new long int[nSpecies];
  nParticleFlag = FALSE;

  theLoads = new oopicList<Load>;
  
  mcc  = _mcc;
  mcti = _mcti;
  dt   = _dt;

  // Set defaults for moving window
  shift = 0;
  shiftDir = 0;
  shiftDelayTime = ( get_grid()->getX()[(95*getJ())/100][0].e1() 
                     - get_grid()->getX()[0][0].e1() ) * iSPEED_OF_LIGHT;
  next_shift_time = -1.;
  shiftNum = 0;
  shift_dt = (  get_grid()->getX()[2][0].e1()
                - get_grid()->getX()[1][0].e1() ) * iSPEED_OF_LIGHT;

#if defined(MPI_DEBUG) && defined(MPI_VERSION)
  dfout << "MPI_RANK = " << MPI_RANK << "; shift_dt = " << shift_dt
        << "\niSPEED_OF_LIGHT = " << iSPEED_OF_LIGHT
        << "\nx2 = " << get_grid()->getX()[2][0].e1()
        << "\nx1 = " << get_grid()->getX()[1][0].e1()
        << "\ndeltaX = " << (  get_grid()->getX()[2][0].e1()
                             - get_grid()->getX()[1][0].e1() )
	<< "\n";
#endif /* defined(MPI_DEBUG) && defined(MPI_VERSION) */

  // cout << "Default shift delay time = " << shiftDelayTime << endl;
  // cout << "Shift delta-t = " << shift_dt << endl;


  oopicListIter<ParticleGroup>	pg(*particleGroupList[0]);
  pg.restart();
  if (!pg.Done()) maxN = pg()->get_maxN();
  else maxN = 500; //	choose a default size

  t = 0.0;  // otherwise simulation time won't be initialized

  // distribute particle group information to boundaries.
  ApplyToList(setSpeciesArray(speciesArray),*boundaryList,Boundary);
  oopicListIter<Species> spiter(*speciesList);
  for(spiter.restart();!spiter.Done();spiter++)
    speciesArray[spiter()->getID()]=spiter();

  // set the ptrNGD and the ptrMapNGDs to zero as initial values
  ptrNGD = 0;
  ptrMapNGDs = 0;

}

SpatialRegion::~SpatialRegion()
{
  oopicListIter<Boundary> bIter(*boundaryList);
  for (bIter.restart(); !bIter.Done();	bIter++)
    bIter.deleteCurrent();
  delete boundaryList;
  boundaryList = NULL;
  for (int i=0; i<nSpecies; i++)
    {
      oopicListIter<ParticleGroup> pgIter(*particleGroupList[i]);
      for (pgIter.restart(); !pgIter.Done();	pgIter++)
        pgIter.deleteCurrent();
    }
  delete[] nParticles;
  delete[] particleGroupList;
  particleGroupList = NULL;
  delete fields;
  delete grid;
  fields = NULL;
  grid = NULL;
  // delete mcc;
  // delete the list of NGD's
  NGDList.deleteAll(); // delete the list elements and data for the NGDs
}

NGD* SpatialRegion::initNGD(const ostring& _gasType, const ostring& _analyticF, 
                            const Scalar& _gasDensity, 
                            const Vector2& _p1Grid, const Vector2& _p2Grid, 
                            const int& discardDumpFileNGDDataFlag)
  throw(Oops){
  // check first if the list is empty:
  if( NGDList.isEmpty() ) {
    try{
      ptrNGD = new NGD(grid, _gasType, _analyticF, _gasDensity,  
                     _p1Grid, _p2Grid, discardDumpFileNGDDataFlag);
    }
    catch(Oops& oops2){
      oops2.prepend("SpatialRegion::initNGD: Error: \n");//ok
      throw oops2;
    }

    NGDList.add(ptrNGD);
    return ptrNGD; 
  }
  
  oopicListIter<NGD> NGDIter(NGDList);
  for ( NGDIter.restart(); !NGDIter.Done(); NGDIter++ ){
    if ( (NGDIter.current())->getGasType() == _gasType ) {
      // this gas density already exists => return a pointer to it.
      return NGDIter.current(); 
    }
  }
  // the GasType was not found in the list, proceed to create it
  try{
    ptrNGD = new NGD(grid, _gasType, _analyticF, _gasDensity,  
                   _p1Grid, _p2Grid, discardDumpFileNGDDataFlag);
  }
  catch(Oops& oops2){
    oops2.prepend("SpatialRegion::initNGD: Error: \n");//OK
    throw oops2;
  }

  NGDList.add(ptrNGD);
  return ptrNGD; 
}

/**
 * A function to give a pointer to the NGDList to the Boundary
 * objects and set the number of NGDs initialized in the 
 * initialized MapNGDs object. 
 */
void SpatialRegion::setNGDListPtrs() {
  // 
  // instantiate an object to handle the NGD mappings 
  // across boundaries.
  // 
  ptrMapNGDs = new MapNGDs(&NGDList, this);
  ptrMapNGDs->setNumNGDs();
  // 
  // allow each Boundary object to point to the MapNGDs object
  // 
  oopicListIter<Boundary>	nextb(*boundaryList);
  for(nextb.restart(); !nextb.Done(); nextb++)
    nextb.current()->setPtrMapNGDs(ptrMapNGDs);
}

#ifdef undef
//--------------------------------------------------------------------
//	Compute the number of particles in this SpatialRegion, return the
// total number, save the number for each species

long SpatialRegion::nParticles(Scalar* nPtcl)
{
  long nTotal = 0;
  for (int i=0; i<nSpecies; i++)
    {
      long n = 0;
      n = CountParticles(i);
      if (nPtcl != NULL) nPtcl[i] = n;
      nTotal += n;
    }
  return nTotal;
}
#endif

long SpatialRegion::CountParticles(int i)
{
  if (nParticleFlag) return nParticles[i]; // caching scheme
  long nTotal = 0;
  long n = 0;
  oopicListIter<ParticleGroup> pgIter(*particleGroupList[i]);
  for (pgIter.restart(); !pgIter.Done(); pgIter++){
    if (pgIter()->get_n()) n += pgIter()->get_n();
    else pgIter.deleteCurrent();
  }
  nTotal += n;
  nParticles[i] = nTotal;
  return nTotal;
}

void SpatialRegion::globalParticleLimit()
{
  long nTotal = 0;
  for (int i=0; i<nSpecies; i++) nTotal += CountParticles(i);
  if (nTotal > particleLimit) increaseParticleWeight();
}

void SpatialRegion::advance() 
  throw(Oops){

#if defined(MPI_DEBUG) && defined(MPI_VERSION)
  dfout << "MPI_RANK = " << MPI_RANK << "; shift_dt = " << shift_dt
	<< endl;
  dfout << "Begin SpatialRegion::advance(): MPI_RANK = " << MPI_RANK 
        << "t = " << t << endl; 
#endif /* defined(MPI_DEBUG) && defined(MPI_VERSION) */

  // Remove the particles that drop below threshold
  if (getBoltzmannFlag()){  
    oopicListIter<Species> siter(*speciesList);
    for (siter.restart(); !siter.Done(); siter++)
      if (siter()->get_collisionModel()==test){
        int i = siter()->getID();
        addParticleList(fields->testParticleList(*particleGroupList[i]));
      }
  }

  // Advance the particles
  try{
    particleAdvance();
  }
  catch(Oops& oops3){
    oops3.prepend("SpatialRegion::advance: Error: \n"); //done
    throw oops3;
  }

  // Advance the fields
  if (fields->get_freezeFields()) t += dt*getFieldSubFlag();
  else fieldsAdvance();
#ifdef MPI_VERSION
  extern MPI_Comm XOOPIC_COMM;
#endif

  // Shift the region
  if ( next_shift_time < 0.0 ) {
    if ( shiftDelayTime <= t ) next_shift_time = t + shift_dt;
    else                       next_shift_time = shiftDelayTime;
  }
  /************
  cout.precision(15);
  cout << "shift_dt        = " << shift_dt << endl
       << "shiftDelayTime  = " << shiftDelayTime << endl
       << "next_shift_time = " << next_shift_time << endl
       << "dt              = " << dt << endl;
  *************/

#if defined(MPI_DEBUG) && defined(MPI_VERSION)
  cout << "MPI_RANK = " << MPI_RANK << "; shift_dt = " << shift_dt << endl;
  cout << "Begin SpatialRegion::advance(): MPI_RANK = " << MPI_RANK 
       << "; t = " << t << endl; 
#endif /* defined(MPI_DEBUG) && defined(MPI_VERSION) */

  if( shiftDir != 0 && (t >= next_shift_time) ) {


#if defined(MPI_DEBUG) && defined(MPI_VERSION)
    cout << "MPI_RANK = " << MPI_RANK << "  t = " << t << "." << endl
         << "next_shift_time = " << next_shift_time << endl
         << "MPI_RANK = " << MPI_RANK << ": shifting region." << endl;
#endif
    try{
      shiftRegion();
    }
    catch(Oops& oops2){
      oops2.prepend("SpatialRegion::advance: Error: \n");//done
      throw oops2;
    }
  }

}

//--------------------------------------------------------------------
//	Advance particles in spatial region in time by dt, including all its
//

void SpatialRegion::particleAdvance()
  throw(Oops){
  nParticleFlag = FALSE;
  fields->clearI(); 
  
  oopicListIter<Species> siter(*speciesList);
  for (siter.restart(); !siter.Done(); siter++){
    Species* species = siter.current();
    int i = species->getID();
    if (species->cycleNow()){			// this species gets pushed now
      if (!fields->get_freezeFields()) {
        fields->clearI(i);
        fields->setAccumulator(i);
      }



      oopicListIter<ParticleGroup> nextpg(*particleGroupList[i]);
      Scalar species_dt = species->get_subcycleIndex()*getFieldSubFlag()*dt;
      species_dt /= species->get_supercycleIndex();
      for (int super=0; super < species->get_supercycleIndex(); super++){
        if (CountParticles(i)>(species->get_particleLimit()))
          increaseParticleWeight(i);
	
  
			
        species->zeroKineticEnergy();
        for (nextpg.restart(); !nextpg.Done(); nextpg++)
          nextpg()->advance(*fields, species_dt);



      // need to pass particles across SRB's
        ApplyToList(passParticles(),*boundaryList,Boundary);
      
        try{
          emit(species); // don't move this without moving the ApplyToList above.
        }
        catch(Oops& oops2){
          oops2.prepend("SpatialRegion::particleAdvance: Error: \n");//done
          throw oops2;
        }
      
        if (mcc) {
          try{
            addParticleList(mcc->collide(*particleGroupList[i]));
          }
          catch(Oops& oops3){
            oops3.prepend("SpatialRegion::particleAdvance: Error: \n");//done
            throw oops3;
          }
        }
      }
      species->resetSubcycleCount();
    }
    else {
      species->incSubcycleCount();
      if (!fields->get_freezeFields()) fields->addtoI(i);
    }

  }


  /**
   * Do the tunneling ionization if the mcti was instantiated. 
   * Note that this call has different logic than the corresponding
   * call for the mcc class given in the species loop above. 
   * The reason is that we need to control the order in which
   * we ionize the different ion species. For this reason we
   * pass directly the pointer to the array of particle group 
   * lists and do the work in derived classes for the different
   * elements for which the tunneling ionization was implemented.