boltzman.cpp

pic 模拟程序！面向对象

  for (j=1;j<=J;j++) 
    for(k=1;k<=K;k++)
      {
	if((bPtr=bcNode[j][k]))
	  {
	    if((bPtr=bc1[j-1][k]))
	      current = Collect(bPtr, .5*dS[j-1][k].e2(), 
				Vector2(j,k),  BoltzmannRho[j][k], scale);
				
	    if((bPtr=bc2[j][k-1]))
	      current += Collect(bPtr, .5*dS[j][k-1].e1(),  
				 Vector2(j,k), BoltzmannRho[j][k], scale);

	    PartCurrent += current;
	    EnergyCurrent += phi[j][k]*current;
	  }
      }

  for (j=1;j<J;j++) 
    for(k=1;k<K;k++)
      {
	if((bPtr=bcNode[j][k]))
	  {
	    if((bPtr=bc1[j][k]))
	      current = Collect(bPtr,.5*dS[j][k].e2(), 
				Vector2(j,k), BoltzmannRho[j][k],scale);
				
	    if((bPtr=bc2[j][k]))
	      current += Collect(bPtr, .5*dS[j][k].e1(),  
				 Vector2(j,k), BoltzmannRho[j][k], scale);
				
	    PartCurrent += current;
	    EnergyCurrent += phi[j][k]*current;
	  }
      }
	
  j=0; k=0;
  if((bPtr=bcNode[j][k]))
    {
      if((bPtr=bc1[j][k]))
	current = Collect(bPtr, .5*dS[j][k].e2(),  
			  Vector2(j,k), BoltzmannRho[j][k], scale);

      if((bPtr=bc2[j][k]))
	current += Collect(bPtr, .5*dS[j][k].e1(),  
			   Vector2(j,k), BoltzmannRho[j][k], scale);

      PartCurrent += current;
      EnergyCurrent += phi[j][k]*current;
    }

  j=0;
  for(k=1;k<K;k++)
    {
      if((bPtr=bcNode[j][k]))
	{
	  if((bPtr=bc1[j][k]))
	    current = Collect(bPtr, .5*dS[j][k].e2(),  
			      Vector2(j,k), BoltzmannRho[j][k], scale);
	  
	  if((bPtr=bc2[j][k]))
	    current += Collect(bPtr, .5*dS[j][k].e1(),  
			       Vector2(j,k), BoltzmannRho[j][k], scale);

	  if((bPtr=bc2[j][k-1]))
	    current += Collect(bPtr, .5*dS[j][k-1].e1(),  
			       Vector2(j,k), BoltzmannRho[j][k], scale);

	  PartCurrent += current;
	  EnergyCurrent += phi[j][k]*current;
	}
    }
	
  if((bPtr=bcNode[j][k]))
    if((bPtr=bc2[j][k-1]))
      {
	current = Collect(bPtr, .5*dS[j][k-1].e1(),  
			  Vector2(j,k), BoltzmannRho[j][k], scale);

	PartCurrent += current;
	EnergyCurrent += phi[j][k]*current;
      }
  
  k=0;
  for (j=1;j<J;j++)
    { 
      if((bPtr=bcNode[j][k]))
	{
	  if((bPtr=bc1[j][k]))
	    current = Collect(bPtr, .5*dS[j][k].e2(),  
			      Vector2(j,k), BoltzmannRho[j][k], scale);

	  if((bPtr=bc1[j-1][k]))
	    current = Collect(bPtr, .5*dS[j-1][k].e2(),  
			      Vector2(j,k), BoltzmannRho[j][k], scale);

	  if((bPtr=bc2[j][k]))
	    current = Collect(bPtr, .5*dS[j][k].e1(),  
			      Vector2(j,k), BoltzmannRho[j][k], scale);
		
	  PartCurrent += current;
	  EnergyCurrent += phi[j][k]*current;
	}
    }

  if((bPtr=bcNode[j][k]))
    if((bPtr=bc1[j-1][k]))
      {
	current = Collect(bPtr, .5*dS[j-1][k].e2(),  
			  Vector2(j,k), BoltzmannRho[j][k], scale);

	PartCurrent += current;
	EnergyCurrent += phi[j][k]*current;
      }
	
  // do correct normalization for particle and energy current

  PartCurrent *= .5*sqrt(2*temperature*fabs(ELECTRON_CHARGE)/M_PI/BoltzSpecies->get_m());  // don't use q_over_m because fabs(ELECTRON_CHARGE) is the converstion from eV to J

  EnergyCurrent *= .5*sqrt(2*temperature/fabs(ELECTRON_CHARGE)/M_PI/BoltzSpecies->get_m());  // don't use q_over_m because fabs(ELECTRON_CHARGE) is the converstion from eV to J

  return (Vector2(PartCurrent, EnergyCurrent)); 	
}

Scalar Boltzmann::Collect(Boundary* bPtr, Scalar area, 
			  const Vector2& x,const Scalar& edge_rho, 
			  const Scalar& scale)
{
  Scalar Flux=0;
  
  Flux = area*edge_rho*bPtr->getBoltzmannFlux();
  if (Flux)
    {
      Scalar periodicFactor=1;
      if (grid->getPeriodicFlagX1()&&(x.e1()==0||x.e1()==J))
	periodicFactor = .5;
      if (grid->getPeriodicFlagX2()&&(x.e2()==0||x.e2()==K))
	periodicFactor = .5;
      
      
      bPtr->collect(*new Particle(x, Vector3(0,0,0), 
				  BoltzSpecies, 
				  fabs(periodicFactor*Flux*scale),1),
		    *new Vector3(0,0,0));  // varweight particle  
    }
  
  return (Flux);
}

ParticleList& Boltzmann::testParticleList(ParticleGroupList& pgList)
{

  if (!pgList.head) return pList;
  Species* species = pgList.head->data->get_species();
  ParticleGroup* pg;
  Scalar phi0;

  getMinPot();

  oopicListIter<ParticleGroup> pgIter(pgList);
  for (pgIter.restart(); !pgIter.Done(); pgIter++)
    {
      pg = pgIter();
      for (int index=0; index<pgIter()->get_n(); index++)
	{
	  //bilinear interpolation
	  phi0 = grid->interpolateBilinear(phi, pg->get_x(index));
	  
	  if ((MinPot-phi0+pg->energy_eV(index))<species->get_threshold())
	    {
	      // make the Boltzmann pg varweight.
	      pList.add(new Particle(pg->get_x(index), pg->get_u(index), 
				     fields->get_bSpecies(), 
				     pg->get_np2c(index), 1));  
				
	      //	Move last particle into this slot
	      pg->fill(index);
	      index--;
	    }
	}
    }
  return pList;
}

void Boltzmann::IncBoltzParticles(ParticleGroupList& pgList)
{

  if (!pgList.head) return;

  Scalar charge = 0;
  Scalar KEint = 0;
  Scalar phi0=0;
  Scalar Denergy = 0;
  Scalar KE = 0;
  
  ParticleGroup* pg;

  int i = BoltzSpecies->getID();
  
  getMinPot();

  oopicListIter<ParticleGroup> pgIter(pgList);
  for (pgIter.restart(); !pgIter.Done(); pgIter++)
    { 
      pg = pgIter();
      for (int index=(pgIter()->get_n()-1); index>=0; index--)
	{
	  phi0 = grid->interpolateBilinear(phi, pg->get_x(index));
	  charge = pg->get_q(index);
	  
	  KE = pg->energy_eV(index);
	  Denergy = charge*(-phi0+pg->energy_eV(index))/qMKS; 
	  Scalar FE = charge*(MinPot-phi0+pg->energy_eV(index))/qMKS; 
	  energy += Denergy;
	  totalCharge += charge;

	}
      pg->set_n(0);
    }  
}

void Boltzmann::getMinPot(void)
{
  // if sign == 1 return a min
  // if sign == -1 return a max
  // need to add something here to see if phi has changed
  //	Scalar MinPot=sign*10000;
  MinPot = q*10000;
  int sign = (int)(q/fabs(q));
  
  for(int j=0;j<=J;j++) for(int k=0;k<=K;k++)
    if (grid->PlasmaRegion(j,k))
      if ((sign*MinPot) > (sign*phi[j][k]))
	MinPot = phi[j][k];	
}


#ifdef HAVE_HDF5

int Boltzmann::dumpH5(dumpHDF5 &dumpObj)
{

  cerr << "entered boltzman::dumpH5(dumpObj)\n";

  hid_t fileId, groupId, datasetId,dataspaceId;
  herr_t status;
  hsize_t rank; 
  hsize_t     phiDim[2];   
  hsize_t     dims;

  hid_t  attrdataspaceId,attributeId;
  hid_t scalarType = dumpObj.getHDF5ScalarType();
  string datasetName;

  // Open the hdf5 file with write access
  fileId = H5Fopen(dumpObj.getDumpFileName().c_str(), 
		   H5F_ACC_RDWR, H5P_DEFAULT);
  groupId = H5Gcreate(fileId,"Boltzmann",0);    

  // dump total charge as attribute
  dims = 1;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);

  attributeId = H5Acreate(groupId, "totalCharge",scalarType, 
			  attrdataspaceId, H5P_DEFAULT);

  status = H5Awrite(attributeId, scalarType, &totalCharge);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);
  
  // dump energy as attribute
  dims = 1;
  attrdataspaceId = H5Screate_simple(1, &dims, NULL);

  attributeId = H5Acreate(groupId, "energy",scalarType, 
			  attrdataspaceId, H5P_DEFAULT);

  status = H5Awrite(attributeId, scalarType, &energy);
  status = H5Aclose(attributeId);
  status = H5Sclose(attrdataspaceId);
  
  datasetName = "phi";
  rank = 2;
  phiDim[0] = J+1;
  phiDim[1] = K+1;
  dataspaceId = H5Screate_simple(rank, phiDim, NULL);

  datasetId = H5Dcreate(groupId, datasetName.c_str(), scalarType, dataspaceId,
			H5P_DEFAULT);
  
  status = H5Dwrite(datasetId, scalarType, H5S_ALL, H5S_ALL,
		    H5P_DEFAULT, phi);
  
  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  H5Gclose(groupId);
  H5Fclose(fileId);
  
  return (1);
}

#endif

int Boltzmann::Dump(FILE* DMPFile)
{

  XGWrite(&totalCharge, ScalarInt, 1, DMPFile, ScalarChar);
  XGWrite(&energy, ScalarInt, 1, DMPFile, ScalarChar);
  for (int j=0; j<=J;j++)
    for (int k=0; k<=K;k++)
      XGWrite(&phi[j][k], ScalarInt, 1, DMPFile, ScalarChar);
  
  return (1);
}


#ifdef HAVE_HDF5

int Boltzmann::restoreH5(dumpHDF5 &restoreObj)
{
  hid_t fileId, groupId, datasetId=0, dataspaceId;
  herr_t status;
  string outFile;
  hid_t scalarType = restoreObj.getHDF5ScalarType();
  hid_t attrId; // attrspaceId,
  
  dump = true;
  
  fileId = H5Fopen(restoreObj.getDumpFileName().c_str(), 
		   H5F_ACC_RDONLY, H5P_DEFAULT);

  groupId = H5Gopen(fileId, "/Boltzmann");
  
  // Read the totalCharge attribute
  attrId = H5Aopen_name(datasetId, "totalCharge");
  status = H5Aread(attrId,scalarType , &totalCharge);
  status = H5Aclose(attrId);
  
  // Read the energy attribute
  attrId = H5Aopen_name(datasetId, "energy");
  status = H5Aread(attrId,scalarType , &energy);
  status = H5Aclose(attrId);
  
  // Read phi  
  datasetId = H5Dopen(groupId,"phi" );
  dataspaceId = H5Dget_space(datasetId );
  
  status = H5Dread(datasetId, scalarType, H5S_ALL, dataspaceId, 
		    H5P_DEFAULT, phi);
  
  H5Sclose(dataspaceId);
  H5Dclose(datasetId);
  H5Gclose(groupId);
  H5Fclose(fileId);

  return (1);
}
#endif

int Boltzmann::Restore(FILE* DMPFile)
{
  dump = true;
  
  XGRead(&totalCharge, ScalarInt, 1, DMPFile, ScalarChar);
  XGRead(&energy, ScalarInt, 1, DMPFile, ScalarChar);
  for (int j=0; j<=J;j++)
    for (int k=0; k<=K;k++)
      XGRead(&phi[j][k], ScalarInt, 1, DMPFile, ScalarChar);
  
  return (1);
}

int Boltzmann::NonlinearSolve(Scalar** source)
{
  int j,k;
  Scalar res=10;
  int it=0;
  int it_out = 0;
  int itmax=20;

  int EB;
  Scalar pstol = .1f;  //tolerance for the poisson solver => tolerance*pstol
  
  getMinPot();  // Finds the min pot within the plasma
  UpdateBRho();
  
  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
    totalRho[j][k] = source[j][k] + BoltzmannRho[j][k];
  
  EBCFlag = 0;
  PSCFlag = 0;
  
  while ((EBCFlag==0)&&(PSCFlag==0)&&(it_out < (itmax/5)))
    {
      while ((EBCFlag==0)&&(PSCFlag==0)&&(it < itmax))
	{	  
	  PSolveCoeff(source); // change psolve coefficients
	  psolve->solve(phi,rho,itermax,pstol*tolerance); 
			
	  getMinPot();
	  UpdateBRho();

	  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
	    totalRho[j][k] = source[j][k] + BoltzmannRho[j][k];
		
	  res = Residual();
	  
	  if (res < tolerance)
	    PSCFlag = 1;

	  EB = EnergyBalance();
	  if (EB==1)
	    {
	      for(j=0;j<=J;j++) for(k=0;k<=K;k++)
		totalRho[j][k] = source[j][k] + BoltzmannRho[j][k];
		
	      res = Residual();
	      it++;
			
	    }
	  
	  else if (EB==2)
	    return (0);
	  
	}
      if (res > tolerance)
	{
	  cout << "Nonlinear Boltzmann Solve FAILURE: Residual = " 
	       << res << endl;
		
	  pstol *=.5;
	  it_out++;
	  cout << "Decreasing psolve tolerance to " << pstol*tolerance << endl;
	  it = 0;
	}
    }
	
  if (res > tolerance)
    {
      cout << "Nonlinear Boltzmann Solve TOTAL FAILURE: Residual = " 
	   << res << endl;
    }
	
  return (1);
}

void Boltzmann::LinearSolve(Scalar** source)
{
	int j,k;
	
	for(j=0;j<=J;j++) for(k=0;k<=K;k++) 
		rho[j][k] = -source[j][k];
	
	switch(ElectrostaticFlag) {
	default:
		fprintf(stderr,"WARNING, poisson solve defaulting to Multigrid in Boltzmann Solve\n");
		
	case ELECTROMAGNETIC_SOLVE:  
	case DADI_SOLVE:
	case PERIODIC_X1_DADI:
		for(j=0;j<=J;j++) for(k=0;k<=K;k++){
			if (b1_orig[j][k])
				b1_local[j][k] = b1_orig[j][k];
			if (b2_orig[j][k])
				b2_local[j][k] = b2_orig[j][k];
		}
		break;
	case CG_SOLVE:
	case GMRES_SOLVE:
		fprintf(stderr,"WARNING, CG and GMRES solves has not been written for Boltzmann Solve\n");
		break;
	case MGRID_SOLVE:
		psolve->PSolveBoltzCoeff(0, qbytemp, phi, 0);
	}

	psolve->solve(phi,rho,itermax,tolerance);

	for(j=0;j<=J;j++) for(k=0;k<=K;k++){
		BoltzmannRho[j][k] = 0.0;
		totalRho[j][k] = source[j][k] + BoltzmannRho[j][k];
	}

}

void Boltzmann::UpdateBRho()
{
	int j,k;
	Scalar sum_of_exp = 0.0;

	for(j=0;j<=J;j++) 
		for(k=0;k<=K;k++){
			if (grid->PlasmaRegion(j,k)){
				if ((-(phi[j][k]-MinPot)*qbytemp)<-50.0)
					BoltzmannRho[j][k] = 0;
				else
					BoltzmannRho[j][k] = exp(-(phi[j][k]-MinPot)*qbytemp);
				sum_of_exp += (BoltzmannRho[j][k]*grid->get_halfCellVolumes()[j][k]);
			}
			else 
				BoltzmannRho[j][k] = 0.0;
		}

	n0 = totalCharge/sum_of_exp;
	
	for(j=0;j<=J;j++) for(k=0;k<=K;k++)
		BoltzmannRho[j][k] *= n0;
	
	return;
}

// int Boltzmann::EnergyBalance(void) -- returns one if balanced
int Boltzmann::EnergyBalance()  
{
  int j,k;

  Scalar PotEnergy = 0.0;
  Scalar PotEnergyDiv = 0.0;
  Scalar PEgrid = 0.0;
  Scalar deltaT = 0.0;

  Scalar OldTemp = temperature;
	
  if (EB_flag)
    {
      for(j=0;j<=J;j++) for(k=0;k<=K;k++)
	{ 
	  PEgrid =  
	    phi[j][k]*BoltzmannRho[j][k]*grid->get_halfCellVolumes()[j][k];
		
	  PotEnergy += PEgrid;
	  PotEnergyDiv += PEgrid*(phi[j][k]-MinPot);
	}
	
      PotEnergy /= fabs(ELECTRON_CHARGE); // J->eV
      PotEnergyDiv *= qbytemp/temperature;
      PotEnergyDiv /= fabs(ELECTRON_CHARGE); // J->eV

      if (PSCFlag&&(fabs((totalCharge*3.0*temperature/qMKS/2.0
			  + PotEnergy - energy)/energy)<tolerance))
	{
	  EBCFlag = 1;
	}
      else
	{
	  PSCFlag = 0;
	  Scalar NewTemp = 2.0*(energy - PotEnergy)*qMKS/totalCharge/3.0;
	  temperature = (NewTemp+temperature)/2;

	  if (temperature<=0)
	    {  
	      // if the temperature goes negative there is not
	      // enought thermal energy for this method to converge
		
	      cout << "temperature went negative.  " << temperature << endl;

	      qbytemp = -1;
	      totalCharge =0;
	      energy = 0;
	      temperature = 0;
	      return 2;

	    }
	  else
	    {
	      qbytemp = q/temperature;
	      EBCFlag = 0;
	      return 0;
	    }
	}
    }
  else
    {
      EBCFlag = 1;
    }
	
  return 1;
}

Vector2 Collisions(Scalar dt)
{	  
  // currently, charge & energy xfer due to collisions
  // not accounted for (but should be?), JH, Sept. 7, 2005
  return Vector2(0,0);
}