fields.h

pic 模拟程序！面向对象

	Vector3** getiL() {return iL;}                 // Added 6-2-95 by KC
	Vector3	getIMesh(int j, int k) {return I[j][k];};
	void	clearI(int i = -1);			//	set Ispecies[i] to 0, -1 means use I
	void	setAccumulator(int i);			// i = species index to accumulate
	void	addtoI(int i);						// for sub_cycling particles
	void	set_iC1(int j, int k, Scalar value) {iC[j][k].set_e1(value);};
	void	set_iC2(int j, int k, Scalar value) {iC[j][k].set_e2(value);};
	void	set_iC3(int j, int k, Scalar value) {iC[j][k].set_e3(value);};
	void	set_iL1(int j, int k, Scalar value) {iL[j][k].set_e1(value);};
   void	set_iL2(int j, int k, Scalar value) {iL[j][k].set_e2(value);};
   void	set_iL3(int j, int k, Scalar value) {iL[j][k].set_e3(value);};
   void	Excite(); //ntg8-16 excites a field component at t = 0
   void	SeedFieldEmitter();
   void  setdt(Scalar dt) {local_dt = dt;};
   BOOL  setBNodeStatic(Vector3 B0,Scalar betwig,Scalar zoff, char* Bfname,const ostring &B01anal,const ostring &B02anal,const ostring &B03anal);   //	static uniform initial magnetic field
   void	setBeamDump(int j1, int j2);	// set up a magnetic beam dump between j1 and j2

   void setEinit(const ostring &E1init, const ostring &E2init, const ostring &E3init);
   void setBinit(const ostring &B1init, const ostring &B2init, const ostring &B3init);
	void  setrho(int j,int k, Scalar value) {rho[j][k]=value;};
	void  setbgrho(int j, int k, Scalar value) {backGroundRho[j][k]+=value;};
   void  setloaddensity(int isp, int j, int k, Scalar &value) {loaddensity[isp][j][k]+=value;};
   void  setShowInitialDensityFlag(int Flag) {ShowInitialDensityFlag=Flag;};
   void  setBGRhoFlag(int Flag) {BGRhoFlag = Flag;};
   Scalar **getrho() {return rho;};
   Scalar getrho(int j, int k) { return rho[j][k]; }; // added for current source, 8/7/02, JeffH
   Scalar **getQ() {return Charge;};
	void  collect_charge_to_grid(ParticleGroup *group,Scalar **_rho);  // collects the charges from particle group 'group' onto the grid rho
	void  compute_rho(ParticleGroupList **list, int nSpecies);  //computes rho, the charge density
	void ZeroCharge(void);
	void  charge_to_charge_density(int nSpecies);  //converts rho from charge on the grid to charge density on the grid
	void addBGRho(void);    // addes DC background rho
	Scalar **getDivDerror() {return DivDerror;}
	void  updateDivDerror(void);   //      Calculate DivD
	void  updatePhi(Scalar **source,Scalar **dest,Scalar t,Scalar dt)throw(Oops); //      Update the potential
	void  initPhi(void) throw(Oops);  //      Initialize our poisson solver
	void  setEfromPhi(void);  //  Make our electric fields consistent with our poisson solution
	Scalar **getphi() { return Phi; };
	void setPhi(Scalar value,int j,int k) { Phi[j][k]=value; }  //set phi at j,k
	void  DivergenceClean(void);  // Do a divergence clean.  call compute_rho and updateDivDerror (in that order) before calling this.
	
	Vector2	getMKS(int j, int k){return grid->getMKS(j, k);}
	Vector2	getMKS(Vector2 x){return grid->getMKS(x);}
	Vector2	getGridCoords(Vector2 xMKS){return grid->getGridCoords(xMKS);}
	Vector3	differentialMove(const Vector2& x0, const Vector3& v, Vector3& u,
		Scalar dt) {return grid->differentialMove(x0, v, u, dt);}
	Grid*	get_grid(){return grid;}
	Scalar getBz0(){return Bz0;}
	Scalar getBx0(){return Bx0;}
	int	getJ() {return grid->getJ();}
	int	getK() {return grid->getK();}
	int   getSub_Cycle_Iter() {return sub_cycle_iter;}
	Scalar   get_dt() { return local_dt; }
	void InodeToI(); //  average the inode currents to the I currents

	void MarderCorrection(Scalar dt);  // Do a marder correction pass

	void setMarderIter(int flag) { MarderIter=flag;}
	void setDivergenceCleanFlag(int flag) { DivergenceCleanFlag=flag;}
	void setCurrentWeightingFlag(int flag) { CurrentWeightingFlag=flag;}
	void setElectrostaticFlag(int flag) { ElectrostaticFlag=flag;}
	void setSynchRadiationFlag(int flag) { SynchRadiationFlag=flag;}
	void set_freezeFields(int flag) {freezeFields=flag;}
	void setNonRelativisticFlag(int flag) { NonRelativisticFlag=flag;}
   void set_nSmoothing(int n) {nSmoothing = n;}
	void setInfiniteBFlag(int flag) { InfiniteBFlag=flag;}
	void setFieldSubFlag(int flag) { FieldSubFlag=flag;}
	void setEMdamping(Scalar value) {EMdamping = value;}
	void setBoltzmannFlag(int flag) { BoltzmannFlag=flag;}
	int getBoltzmannFlag() {return BoltzmannFlag;}
	void setMarderParameter(Scalar value) { MarderParameter = value;}
	int getElectrostatic() { return ElectrostaticFlag; }
	int getSynchRadiation() { return SynchRadiationFlag; }
	int get_freezeFields() {return freezeFields;}
	Scalar getEMdamping() {return EMdamping;}
	int getNonRelativistic() { return NonRelativisticFlag; }
	int getInfiniteBFlag() { return InfiniteBFlag; }
	int getFieldSub() { return FieldSubFlag; }
	PoissonSolve *getPoissonSolve() { return psolve; }
	void set_epsi(Scalar value,int j,int k) {epsi[j][k]=value;};
	Scalar get_epsi(int j, int k) {return epsi[j][k];};
	void setBoltzmann(Boltzmann *B) { theBoltzmann=B; };
	Scalar *** getRhoSpecies() { return rho_species; };
	Scalar ** getBoltzmannRho() {return theBoltzmann->getBoltzmannRho();};
	Species* get_bSpecies() {return bSpecies;};
	void set_bSpecies(Species *BS) {bSpecies=BS;};
	void CopyBoltzmannRho() {theBoltzmann->CopyBoltzmannRho();}
	Scalar ** getTotalRho() {return theBoltzmann->getTotalRho();};
	Scalar ** getdelPhi() { return delPhi; };
	Scalar getPhiP(int j, int k) { return PhiP[j][k]; }; // added for current source, 8/7/02, JeffH
	Scalar *** getloaddensity() { return loaddensity;};
	int getShowInitialDensityFlag() {return ShowInitialDensityFlag;}; 
	void compute_iC_iL();
//  Scalar getMinPot(int sign);
//	void setBoltzQE(Scalar charge, Scalar energy) {theBoltzmann->setQE(charge,energy);};
	ParticleList& testParticleList(ParticleGroupList& pgList) {return theBoltzmann->testParticleList(pgList);};
	void IncBoltzParticles(ParticleGroupList& pgList) {theBoltzmann->IncBoltzParticles(pgList);};


#ifdef HAVE_HDF5

	int dumpH5(dumpHDF5 &dumpObj);
	int restoreH5(dumpHDF5 &restoreObj);
#endif

	int Dump(FILE *DMPFile);
	int Restore(FILE *DMPFile);
	int Restore_2_00(FILE *DMPFile);

  // morgan
  void set_up_inverter(Grid*) throw(Oops);

// support for moving window
// direction:  1 = j+1 ==> j, 2 = j-1 ==> j, 3 = k+1 --> k, 4 = k-1 --> k

/**
 * Send fields in a given direction, which corresponds to a certain
 * boundary
 *
 * @param direction The direction of the shift
 * @param bPtr the boundary that fields are going out
 */
  void sendFieldsAndShift(int direction, Boundary* bPtr);

/**
 * Receive fields due to a shift of a given direction.
 *
 * @param direction The direction of the shift
 * @param bPtr the boundary that fields are going in from
 */
  void recvShiftedFields(int direction, Boundary* bPtr,Scalar t, Scalar dt);  

/**
 * Load in any new fields for moving boundary
 *
 */
  void loadShiftedFields(Scalar t, Scalar dt);  

};

//--------------------------------------------------------------------
//	Accumulate current based on x(t) and x(t+dt*f).  Here frac_v3dt
//	is the fraction of the timestep for this piece of the move times
//	the azimuthal displacement v3*dt.  Cell volume is more accurately
//	chosen separately for each current element.
// q_dt == q/dt

inline void Fields::accumulate(Vector2& xOld, Vector2&x, Scalar frac_v3dt,
	Scalar q_dt)
{
  if (get_freezeFields()) return; // emit always does EM push
	int	j = (int) MIN(xOld.e1(),x.e1());
   int	k = (int) MIN(xOld.e2(),x.e2());
	if (j == grid->getJ()) j--;
	if (k == grid->getK()) k--;
	Vector2	dw = x - xOld;
	Vector2	avg_w = 0.5*(xOld + x) - Vector2(j, k);

	frac_v3dt /= grid->cellVolume(j, k);
	if(CurrentWeightingFlag==0)
	{
          if(accumulator) {
		accumulator[j][k] += q_dt*Vector3(dw.e1()*(1 - avg_w.e2()),
									dw.e2()*(1 - avg_w.e1()),
									grid->dS3Prime(j, k)*(1 - avg_w.e1())*
									(1 - avg_w.e2())*frac_v3dt);

		accumulator[j][k+1] += q_dt*Vector3(dw.e1()*avg_w.e2(),
									  0,
									  grid->dS3Prime(j, k+1)*
									  (1 - avg_w.e1())*avg_w.e2()*frac_v3dt);
		accumulator[j+1][k] += q_dt*Vector3(0,
									  dw.e2()*avg_w.e1(),
									  grid->dS3Prime(j+1, k)*
									  avg_w.e1()*(1 - avg_w.e2())*frac_v3dt);

		accumulator[j+1][k+1] += q_dt*Vector3(0,
										 0,
										 grid->dS3Prime(j+1, k+1)*avg_w.e1()*avg_w.e2()
										 *frac_v3dt);
          }
	}
	else  /*  the bilinear weighting of current */
	{
          if(accumulator) {
		accumulator[j][k] += q_dt*Vector3(dw.e1()*(1 - avg_w.e2())*(1-avg_w.e1()),
										 dw.e2()*(1 - avg_w.e1())*(1-avg_w.e2()),
										 grid->dS3Prime(j, k)*(1 - avg_w.e1())*
										 (1 - avg_w.e2())*frac_v3dt);

		accumulator[j][k+1] += q_dt*Vector3(dw.e1()*avg_w.e2()*(1-avg_w.e1()),
									  dw.e2()*avg_w.e2()*(1-avg_w.e1()),
									  grid->dS3Prime(j, k+1)*
									  (1 - avg_w.e1())*avg_w.e2()*frac_v3dt);

		accumulator[j+1][k] += q_dt*Vector3(dw.e1()*avg_w.e1()*(1-avg_w.e2()),
									  dw.e2()*avg_w.e1()*(1-avg_w.e2()),
									  grid->dS3Prime(j+1, k)*
									  avg_w.e1()*(1 - avg_w.e2())*frac_v3dt);
		accumulator[j+1][k+1] += q_dt*Vector3(dw.e1()*avg_w.e1()*avg_w.e2(),
										 dw.e2()*avg_w.e1()*avg_w.e2(),
										 grid->dS3Prime(j+1, k+1)*avg_w.e1()*avg_w.e2()
										 *frac_v3dt);
          }
	}
}

#endif // #ifndef __FIELDS_H