fields.cpp

pic 模拟程序！面向对象

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

   fields.cpp

   Purpose:	Member functions for Fields class.  This class provides 
   the full set of electromagnetic fields, as well as 
   methods for advancing the PDEs in time.


   Revision/Programmer/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	(JohnV 08-01-94) Add prototypes for dadi() and init_dadi().
   0.980	(JohnV 08-17-95) Add compute_iC_iL(), generalize for arbitrary meshed epsilon.
   0.99  (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.991	(JohnV 09-10-95) Add setBeamDump() to create a magnetic beam dump;
   required making BNodeStatic a Vector3.
   0.992	(Ron Cohen 09-14-95) More parameters for BNodeStatic added to
   allow for more general initial magnetic fields.
   1.001	(JohnV, KC 04-19-96) Optimized innermost loop of advanceE() +3% speed.
   1.002	(JohnV, 04-25-96) Fixed nomenclature in compute_iC_iL().
   1.003 (JohnV 05-03-96) Modified toNodes() to set interior fields only, and call 
   Boundary::toNodes() for edges and internal boundaries.
   1.0035 (KC 5-15-96) Added EMdamping.  Naoki will tells us what we are doing to the 
   physics.
   1.004	(JohnV 05-15-96) Add speciesList, Ispecies to improve subcycling.
   1.01  (JohnV 09-30-97) Add freezeFields, which uses initial fields.
   1.02  (KC 01-13-98) Changed Compute_ic_iL to use dl, dlPrime, dS, and dSPrime from Grid.
   2.01  (JohnV 03-19-98) Add nSmoothing to apply binomial filter to rho within computeRho.
   2.02  (JohnV 02-04-99) MarderFlag->MarderIter 
   2.03  (Bruhwiler 11-08-99) added more vector and scalar quantities 
   to the ShiftFields() method, so that moving window algorithm 
   doesn't break emdamping algorithm
   2.0?  (Cary 22 Jan 00) Broke up shift fields into two parts in anticipation
   of MPI installation.
   CVS1.80.2.6 (Cary 23 Jan 00) Reindented toNodes so that is easily seen on
   80 char telnet window.  No significant changes here.  For 
   sendFieldsAndShift, moved the send to after the shift.  Now 
   detect if boundary is SPATIAL_REGION_BOUNDARY, if so, call
   putCharge_and_Current.  Similarly, in recvShiftedFields, I check
   to see if boundary is SPATIAL_REGION_BOUNDARY.  If so, I call
   applyFields.
   CVS1.80.2.9 (Cary 26 Jan 00) changed sendFieldsAndShift to shiftFieldsAndSend
   for accuracy.  
   CVS1.87 (JohnV 05 Jan 2001) Added check for magnetic field file in setBNodeStatic(),
   improved efficiency.
   CVS1.88 (JohnV 21 Mar 2001) Re-ordered field components to B1, B2, B3 when
   magnetic field file

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

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

#ifdef PETSC
#include "parpoi.h"
#endif

#include "fields.h"
#include "ptclgrp.h"
#include "dadirz.h"
#include "dadixy.h"
#include "dadixyp.h"
#include "boltzman.h"
#include "multigrid.h"
#include "dump.h"

//#include "cgxy.h"
// morgan
#include "domain.h"
#include "electrostatic_operator.h"
#include "cg.h"
#include "gmres.h"

#include "eval.h"
#include <cstdio>

#ifdef _MSC_VER
#include <iomanip>
#include <iostream>
#include <fstream>
using std::cout;
using std::endl;
#else
#include <txc_streams.h>
using namespace std;
#endif

// io includes

#ifdef HAVE_HDF5
#include <hdf5.h>
#include "dumpHDF5.h"
#endif



#ifdef PHI_TEST
Scalar analytic_rho(int j, int k);
#endif
extern Evaluator *adv_eval;


//--------------------------------------------------------------------
//	Allocate arrays for the electromagnetic fields and the capacitance
//	and inductance matrices.

Fields::Fields(Grid* _grid, BoundaryList* _boundaryList, SpeciesList* _speciesList, 
               Scalar epsR, int _ESflag,Scalar _presidue, int _BoltzmannFlag, int _SRflag)
{
   int i, j, k;
   grid = _grid;
   boundaryList = _boundaryList;
   speciesList = _speciesList;
   nSpecies = speciesList->nItems();
   presidue = _presidue;
   // Initialize flags:
   MarderIter=0;
   BGRhoFlag=0;
   DivergenceCleanFlag=0;
   CurrentWeightingFlag=0;
   ElectrostaticFlag=_ESflag;
   SynchRadiationFlag = _SRflag;
   freezeFields = 0;
   BoltzmannFlag = _BoltzmannFlag;
   ShowInitialDensityFlag = FALSE;
   MarderParameter=0;
   psolve=0;
   sub_cycle_iter = 0;
   EMdamping = 0;
   epsilonR = epsR;
   J = grid->getJ();
   K = grid->getK();
   //grid->setBoundaryMask(*boundaryList);
   //	allocate memory
   rho = new Scalar *[J + 1];
   Charge = new Scalar *[J + 1];
   backGroundRho = new Scalar *[J + 1];
   Phi = new Scalar *[J + 1];
   PhiP = new Scalar *[J + 1];
   DivDerror = new Scalar *[J + 1];
   ENode = new Vector3* [J + 1];
   BNode = new Vector3* [J + 1];
   BNodeStatic = new Vector3* [J+1];
   BNodeDynamic = new Vector3* [J+1];
   I = new Vector3* [J + 1];
   Ispecies = new Vector3** [nSpecies];
   oopicListIter<Species> spIter(*speciesList);
   for (spIter.restart(); !spIter.Done(); spIter++){
      i = spIter.current()->getID();
      if (spIter.current()->isSubcycled()) Ispecies[i] = new Vector3* [J+1];
      else Ispecies[i] = NULL;
   }
   Inode = new Vector3* [J + 1];
   intEdl = new Vector3* [J + 1];
   intEdlPrime = new Vector3* [J + 1];
   intEdlBar = new Vector3* [J + 1];
   intBdS = new Vector3* [J + 1];
   iC = new Vector3* [J + 1];
   iL = new Vector3* [J + 1];
   epsi = new Scalar *[J];  //cell centered
   rho_species = 0;
   accumulator = 0;

   for (j=0; j<=J; j++)
      {
         k = 0;
         rho[j] = new Scalar[K + 1];
         Charge[j] = new Scalar[K + 1];
         backGroundRho[j] = new Scalar[K + 1];
         Phi[j] = new Scalar [K + 1];
         PhiP[j] = new Scalar [K + 1];
         DivDerror[j]= new Scalar[K + 1];
         ENode[j] = new Vector3[K + 1];
         BNode[j] = new Vector3[K + 1];
         BNodeStatic[j] = new Vector3[K+1];
         BNodeDynamic[j] = new Vector3[K+1];
         I[j] = new Vector3[K + 1];
         for (i = 0; i < nSpecies; i++)
            if (Ispecies[i]) Ispecies[i][j] = new Vector3[K + 1];
         Inode[j] = new Vector3[K + 1];
         intEdl[j] = new Vector3[K + 1];
         intEdlPrime[j] = new Vector3[K + 1];
         intEdlBar[j] = new Vector3[K + 1];
         intBdS[j] = new Vector3[K + 1];
         iC[j] = new Vector3[K + 1];
         iL[j] = new Vector3[K + 1];
      }
   
   for(j=0; j<J; j++)
      epsi[j] = new Scalar[K];

   loaddensity = new Scalar **[nSpecies];
   for(i=0;i<nSpecies;i++) {
      loaddensity[i] = new Scalar *[J+1];
      for(j=0;j<=J;j++) {
         loaddensity[i][j]=new Scalar[K+1];
         //zero the memory
         memset(loaddensity[i][j],0,(K+1)*sizeof(Scalar));
      }
   }

   intEdlBasePtr = intEdl;
   //  The following initializations to zero are important
   //  correct answers depend upon them.
   for(j=0;j<=J;j++)
      for(k=0;k<=K;k++) {
         rho[j][k]=0;Phi[j][k]=0;backGroundRho[j][k]=0;
         Charge[j][k]=0;
         PhiP[j][k]=0;
         DivDerror[j][k]=0;
         //			epsi[j][k]=1/iEPS0;
      }; //initialize rho to zero
   for(j=0;j<J;j++)
      for(k=0;k<K;k++)
         epsi[j][k]=1/iEPS0;
   minusrho=0;
   d=0; delPhi=0;  //These initializations are important
   oopicListIter<Boundary>	nextb(*boundaryList);
   for(nextb.restart(); !nextb.Done(); nextb++)
      nextb.current()->setFields(*this);
   compute_iC_iL();
   // must be done AFTER setFields()
   initPhi();  // initialize the Poisson solver  MUST be done before the boundaries
   // are done below
   for (nextb.restart(); !nextb.Done(); nextb++)
      nextb.current()->setPassives();
   //  Initialize the Laplace solutions for the potential
   //  due to time-varying equipotential boundaries
#ifdef BCT_DEBUG
   printf("\n");
   for(k=K;k>=0;k--) {
      for(j=0;j<=J;j++) {
         Boundary *B = grid->GetNodeBoundary()[j][k];
         int type;
         if(B!=NULL) type = B->getBCType();
         else type = FREESPACE;

         printf("%d\t",type);
      }
      printf("\n");
   }
#endif
#undef RANDOM_FIELDS
#ifdef RANDOM_FIELDS
   // Here we want to set ALL the fields in the system to something random and predictable.
   for(j=0;j<=J;j++) 
      for(k=0;k<=K;k++) {
         int j1,k1, seed;
         k1 = k;
         j1 = (int)(( (grid->getMKS(j,k).e1())/grid->dl1(1,1)) + 0.5);
         seed = k + j1 * K;
         srand(seed);
         /*			intEdl[j][k] = Vector3((rand()/(float)RAND_MAX - 0.5),(rand()/(float)RAND_MAX - 0.5),(rand()/(float)RAND_MAX - 0.5));
                                intBdS[j][k] = Vector3((rand()/(float)RAND_MAX - 0.5),(rand()/(float)RAND_MAX - 0.5),(rand()/(float)RAND_MAX - 0.5)); */

         intEdl[j][k] = Vector3(seed, 10*seed, 100*seed);
         intBdS[j][k] = 1e-7 * Vector3(.001*seed,.01*seed, .1*seed);
      }
   
#endif			
   

   
   for (nextb.restart(); !nextb.Done(); nextb++)
      nextb.current()->InitializeLaplaceSolution();
}

//--------------------------------------------------------------------
//	Deallocate memory for Fields object

Fields::~Fields()
{
   int i;
   int j;
   for (j=0; j<=J; j++)
      {
         delete[] ENode[j];
         delete[] BNode[j];
         delete[] BNodeStatic[j];
         delete[] BNodeDynamic[j];
         delete[] I[j];
         delete[] Inode[j];
         for (i=0; i < nSpecies; i++) if (Ispecies[i]) delete[] Ispecies[i][j];
         delete[] intEdl[j];
         delete[] intEdlPrime[j];
         delete[] intEdlBar[j];
         delete[] intBdS[j];
         delete[] iC[j];
         delete[] iL[j];
         delete[] DivDerror[j];
         delete[] rho[j];
         delete[] Charge[j];
         delete[] backGroundRho[j];
         delete[] Phi[j];
         delete[] PhiP[j];
         if(delPhi!=0) delete[] delPhi[j] ;
         if(minusrho!=0) delete[] minusrho[j]; 
         if(d!=0) delete[] d[j];
      }
   for (j=0; j<J; j++)
      delete[] epsi[j];
   delete[] ENode;
   delete[] BNode;
   delete[] BNodeStatic;
   delete[] BNodeDynamic;
   delete[] I;
   delete[] Inode;
   for (i=0; i < nSpecies; i++) {
      if(Ispecies[i])
         delete[] Ispecies[i];
      else
         delete Ispecies[i];
   }
   delete[] Ispecies;
   delete[] intEdl;
   delete[] intEdlPrime;
   delete[] intEdlBar;
   delete[] intBdS;
   delete[] iC;
   delete[] iL;
   delete[] DivDerror;
   delete[] rho;
   delete[] Charge;
   delete[] backGroundRho;
   delete[] Phi;
   delete[] PhiP;
   delete[] epsi;
   if(delPhi!=0) delete[] delPhi;
   if(minusrho!=0)   delete[] minusrho;
   if(d!=0) delete[] d;
   delete psolve;
   if(rho_species)
      {
         int J=grid->getJ();
         for(i=0; i<nSpecies; i++)
            {
               for(j=0; j<J+1; j++) 
                  delete [] rho_species[i][j];
               delete [] rho_species[i];
            }
         delete [] rho_species;
      }
   int J=grid->getJ();
   for(i=0; i<nSpecies; i++)
      {
         for(j=0; j<=J; j++) 
            delete [] loaddensity[i][j];
         delete [] loaddensity[i];
      }
   delete [] loaddensity;
   /* begin 
      if (poisson) {
      delete poisson;
      poisson = 0;
      } 
      */
   /* 
      if (dom) {
      delete dom;
      dom = 0;
      }
      
      end debug */
}

//--------------------------------------------------------------------
//	Advance the fields from t-dt->t by solving Maxwell's equations.

void Fields::advance(Scalar t, Scalar dt)
  throw(Oops){

   /*
      int kout = K/4;
      #ifdef MPI_VERSION
      extern int MPI_RANK;
      if(MPI_RANK == 0){
      cout << "MPI_RANK = " << MPI_RANK << ", Fields::advance entered.  kout = "
      << kout << "." << endl;
      #else
      cout << "Fields::advance entered.  kout = " <<
      kout << "." << endl;
      #endif
      cout << "E1 =";
      for(int j1=0; j1<=J; j1++)  cerr << " " << intEdl[j1][kout].e1();
      cout << endl;
      cout << "E2 =";
      for(int j1=0; j1<=J; j1++)  cerr << " " << intEdl[j1][kout].e2();
      cout << endl;
      cout << "I1 =";