boltzman.cpp

pic 模拟程序！面向对象

#include "boltzman.h"
#include "psolve.h"
#include "ovector.h"
#include "misc.h"
#include "ptclgrp.h"
#include "particle.h"
#include "fields.h"
#include <math.h>

#ifdef _MSC_VER
using std::cout;
using std::cerr;
using std::endl;
#endif

/* 
   Will only work with Multigrid and maybe with DADI - 
   Keith will be happy to fix if you have a problem -
   send complaints to [e-mail]
*/

Boltzmann::Boltzmann(Fields *F, Scalar _temperature, 
		     Species* _BoltzSpecies, bool _EB_flag)
{
  dump = false;
  fields = F;
  BoltzSpecies = _BoltzSpecies;
  grid = fields->get_grid();
  psolve = fields->getPoissonSolve();
  ElectrostaticFlag = fields->getElectrostatic();
  bcNode = grid->GetNodeBoundary();

  itermax=100;  //maximum # of iterations for psolve

  //tolerance for the boltzmann solve is set from input file
  tolerance = fields->getpresidue(); 
  
  qMKS = BoltzSpecies->get_q();
  q = qMKS/fabs(ELECTRON_CHARGE);
  J = fields->getJ();
  K = fields->getK();

  temperature = _temperature;
  qbytemp = q/temperature;

  EB_flag = _EB_flag;  // if EB_flag the energy balance is used 
                       // otherwise the temperature is constant

  // allocate field arrays
  foo = new Scalar *[J+1]; 
  phi = new Scalar *[J+1]; 
  rho = new Scalar *[J+1];
  BoltzmannRho = new Scalar *[J+1];
  totalRho = new Scalar *[J+1];
  int j, k;
  for(j=0;j<=J;j++)
    {
      foo[j] = new Scalar[K+1];
      phi[j] = new Scalar[K+1];
      rho[j] = new Scalar[K+1];
      BoltzmannRho[j] = new Scalar[K+1];
      totalRho[j] = new Scalar[K+1];
    }
  
  // zero field arrays
  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
    {
      phi[j][k] = 0.0;
      foo[j][k] = 0.0;
      BoltzmannRho[j][k] = 0.0;
      totalRho[j][k] = 0.0;
    }
  
  switch(ElectrostaticFlag) 
    {
    default:
      fprintf(stderr,
	      "WARNING, poisson solve defaulting in Boltzmann Solve\n");
      
    case ELECTROMAGNETIC_SOLVE:  
    case DADI_SOLVE:
    case PERIODIC_X1_DADI:
      b1_orig = new Scalar *[J+1];
      b2_orig = new Scalar *[J+1]; 
      
      for(j=0;j<=J;j++)
	{
	  b1_orig[j] = new Scalar[K+1];
	  b2_orig[j] = new Scalar[K+1];
	}
      
      // *1_local, *2_local is shared by psolve and boltzmann
      a1_local = psolve->get_a1_coefficients();
      a2_local = psolve->get_a2_coefficients();
      b1_local = psolve->get_b1_coefficients();
      b2_local = psolve->get_b2_coefficients();
      c1_local = psolve->get_c1_coefficients();
      c2_local = psolve->get_c2_coefficients();
      
      for(j=0;j<=J;j++) for(k=0;k<=K;k++)
	{
	  b1_orig[j][k] = b1_local[j][k];
	  b2_orig[j][k] = b2_local[j][k];
	}
      break;
      
    case CG_SOLVE:
      fprintf(stderr,
	      "WARNING, CG solve has not been written for Boltzmann Solve\n");
      exit(1);
      break;
      
    case GMRES_SOLVE:
      fprintf(stderr,
	      "WARNING, GMRES has not been written for Boltzmann Solve\n");
      exit(1);
      break;
      
    case MGRID_SOLVE:
	    
      resCoeff = new Scalar**[J+1];
      for (j=0;j<=J;j++)
	{
	  resCoeff[j] = new Scalar*[K+1];
	  for (k=0;k<=K;k++)
	    resCoeff[j][k] = new Scalar[5];
	}		
      resCoeff_local = psolve->get_resCoeff();

      for(j=0;j<=J;j++) for(k=0;k<=K;k++)
	for (int s=0; s<5; s++)
	  resCoeff[j][k][s] = resCoeff_local[j][k][s];
      break;
    }
}

Boltzmann::~Boltzmann()
{

  // deallocate grid arrays
  for(int j=0; j<=J; j++)
    {
      switch(ElectrostaticFlag) 
	{
	default:
	  fprintf(stderr,
		  "WARNING, poisson solve defaulting in Boltzmann Solve\n");
	  
	case ELECTROMAGNETIC_SOLVE:  
	case DADI_SOLVE:
	case PERIODIC_X1_DADI:	
	  
	  delete[] b1_orig[j];
	  delete[] b1_local[j];	
	  break;
	  
	case MGRID_SOLVE:
	  for(int k=0; k < K; k++)
	    {
	      delete[] resCoeff[j][k];
	    }
	  delete[] resCoeff[j];
	  break;
	}
      delete[] phi[j];
      delete[] BoltzmannRho[j];
      delete[] rho[j];
      delete[] totalRho[j];
      delete[] foo[j];
    }
  
  switch(ElectrostaticFlag) 
    {
    default:
      fprintf(stderr,
	      "WARNING, poisson solve defaulting in Boltzmann Solve\n");
      
    case ELECTROMAGNETIC_SOLVE:  
    case DADI_SOLVE:
    case PERIODIC_X1_DADI:
      
      delete[] b1_orig;
      delete[] b1_local;
      break;
    
    case MGRID_SOLVE:
      delete[] resCoeff;
      break;
    }
  delete[] phi;
  delete[] BoltzmannRho;
  delete[] rho;
  delete[] totalRho;
  delete[] foo;
}

// called from Fields::ElectrostaticSolve()
void Boltzmann::updatePhi(Scalar **source, Scalar **dest, 
			  const Scalar t, Scalar dt)
{
  int j,k;
  static int init_flag=1;
  
  if(EB_flag)
    {
      fprintf(stderr, "Energy Balance flag is on (buggy!)[press enter]\n");
      getchar();
    }

  if (init_flag)
    {
      init(source, dest);
      init_flag = 0;
      return;
    }
  
  // copy Fields phi (dest) to local phi
  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
    phi[j][k] = dest[j][k];
	
  if(totalCharge)
    {
      if (totalCharge*BoltzSpecies->get_q()<0)
	{
	  totalCharge = 0;
	  energy = 0;
	}
    }
  
  if(totalCharge)
    {
      // why is this done?, JH, Aug. 2005
      if (temperature<.05)
	{
	  totalCharge = 0.0;
	  temperature = 6; 
	  energy = 0;
	}
    }      

  if(totalCharge)
    {
      if(NonlinearSolve(source))
	{
	  BoltzSpecies->zeroKineticEnergy();
	  BoltzSpecies->addKineticEnergy(3.0/2.0*temperature*totalCharge/q);
	  
	  // deltaC: change in total charge (.e1())
	  //  and energy (.e2())
	  Vector2 deltaC = Vector2(0,0); 

	  // currently, charge and energy xfer due to collisions
	  // not accounted for (but should be?), JH, Sept. 7, 2005
	  //			deltaC = Collisions(dt);

	  // particle and energy flux to wall
	  deltaC += dt*ParticleFlux(dt);
				  
	  // update total charge and energy quantities
	  totalCharge += deltaC.e1();
	  energy += 5.0/2.0*temperature*deltaC.e1()/qMKS;
	  energy += deltaC.e2();	  
	}
      else
	{
	  LinearSolve(source);
	  BoltzSpecies->zeroKineticEnergy();
	  BoltzSpecies->addKineticEnergy((3.0/2.0*totalCharge*temperature)/q);
	}
    }
  else
    {
      LinearSolve(source);
      BoltzSpecies->zeroKineticEnergy();
    }

  // copy local phi back to Fields:: phi (dest)
  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
    dest[j][k]=phi[j][k];  
}

// initial Boltzmann solver first time updatePhi is called
void Boltzmann::init(Scalar **source, Scalar **dest)
{
  int j,k;
  Scalar PotEnergy = 0;
  Scalar totalVolume = 0;
  bool EBFlag_orig;

  // copy Fields:: phi (dest) to local phi
  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
    if (!phi[j][k])
      phi[j][k] = dest[j][k];
  
  if (!dump)
    {
      // set the init totalCharge equal - particle density
      totalCharge=0.0; 
      for(j=0;j<=J;j++) 
	for(k=0;k<=K;k++) 
	  totalCharge+=source[j][k]*grid->get_halfCellVolumes()[j][k];
      totalCharge *= q; 
      
      LinearSolve(source);
	
      // calculate initial energy
      energy = temperature*3.0/2.0*totalCharge; 

      // swap energy balance flag and zero
      EBFlag_orig = EB_flag;
      EB_flag = false;
      
      // Start off with this temperture and find the PotEnergy
      bRho = fields->getRhoSpecies()[BoltzSpecies->getID()];
      
      if (totalCharge)
	NonlinearSolve(source);
      else
	LinearSolve(source);

      EB_flag = EBFlag_orig; // swap energy balance flag back
    }
  else // if dumped
    {
      bRho = fields->getRhoSpecies()[BoltzSpecies->getID()];
      
      if (totalCharge)
	NonlinearSolve(source);
      else
	LinearSolve(source);
    }

  // copy local phi back to Fields:: phi (dest)
  for(j=0;j<=J;j++) for(k=0;k<=K;k++)
    dest[j][k]=phi[j][k];
}

void Boltzmann::PSolveCoeff(Scalar **source)
{
  
  Scalar temp;
  int j,k;
  
  switch(ElectrostaticFlag) 
    {
    default:
      fprintf(stderr,
	      "WARNING, poisson solve defaulting 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++)
	    {
	      temp = -BoltzmannRho[j][k];
	      //		temp=-n0*exp(-phi[j][k]*qbytemp);
	      rho[j][k] = -source[j][k]+(1+qbytemp*phi[j][k])*temp;
	      if (b1_orig[j][k])
		b1_local[j][k] = b1_orig[j][k]+0.5*qbytemp*temp;
	      if (b2_orig[j][k])
		b2_local[j][k] = b2_orig[j][k]+0.5*qbytemp*temp;
	    }
	  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(n0, qbytemp, BoltzmannRho, MinPot);
      for(j=0;j<=J;j++) for(k=0;k<=K;k++){
	temp = -BoltzmannRho[j][k];
	rho[j][k] = -source[j][k]+(1+qbytemp*phi[j][k])*temp;
      }
      break;
    }
}

Scalar Boltzmann::Residual()
{
  Scalar res = 0.0;
  Scalar temp = 0.0;
	
  switch(ElectrostaticFlag) 
    {
    default:
      fprintf(stderr,
	      "WARNING, poisson solve defaulting in Boltzmann Solve\n");
		
    case ELECTROMAGNETIC_SOLVE:  
    case DADI_SOLVE:
      int j,k;
		
      j=0;k=0;
      if (b1_orig[j][k] !=0)  
	{
	  temp += totalRho[j][k];
	  res += sqr(b1_orig[j][k]*phi[j][k]+c1_local[j][k]*phi[j+1][k]+
		     b2_orig[j][k]*phi[j][k]+c2_local[j][k]*phi[j][k+1]+
		     totalRho[j][k]);
	  
	  foo[j][k] = sqr(b1_orig[j][k]*phi[j][k]+c1_local[j][k]*phi[j+1][k]+
			  b2_orig[j][k]*phi[j][k]+c2_local[j][k]*phi[j][k+1]+
			  totalRho[j][k]); 
	}

      k=K;
      if (b1_orig[j][k] !=0)
	{
	  temp += totalRho[j][k];
	  res += sqr(b1_orig[j][k]*phi[j][k]+c1_local[j][k]*phi[j+1][k]+
		     a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]+
		     totalRho[j][k]);
	  
	  foo[j][k] = sqr(b1_orig[j][k]*phi[j][k]+c1_local[j][k]*phi[j+1][k]+
			  a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]+
			  totalRho[j][k]);  
	}

      j=J;
      if (b1_orig[j][k] !=0)
	{
	  temp += totalRho[j][k];
	  res += sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]+
		     a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]+
		     totalRho[j][k]);
	  
	  foo[j][k] = sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]+
			  a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]+
			  totalRho[j][k]); 
	}
		
      k=0;
      
      if (b1_orig[j][k] !=0)
	{
	  temp += totalRho[j][k];
	  res += sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]+
		     b2_orig[j][k]*phi[j][k]+c2_local[j][k]*phi[j][k+1]+
		     totalRho[j][k]);
	  
	  foo[j][k] = sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]+
			  b2_orig[j][k]*phi[j][k]+c2_local[j][k]*phi[j][k+1]+
			  totalRho[j][k]);
	}

      for (j=1;j<J;j++)
	{
	  k = 0;
	  if (b1_orig[j][k] !=0){
	    temp += totalRho[j][k];
	    res += sqr(a1_local[j][k]*phi[j-1][k]+
		       b1_orig[j][k]*phi[j][k]+c1_local[j][k]*phi[j+1][k]+
		       b2_orig[j][k]*phi[j][k]+c2_local[j][k]*phi[j][k+1]+
		       totalRho[j][k]);

	    foo[j][k] = sqr(a1_local[j][k]*phi[j-1][k]
			    +b1_orig[j][k]*phi[j][k]
			    +c1_local[j][k]*phi[j+1][k]+
			    b2_orig[j][k]*phi[j][k]
			    +c2_local[j][k]*phi[j][k+1]+
			    totalRho[j][k]);
	  }
	  
	  k = K;
	  if (b1_orig[j][k] !=0)
	    {
	      temp += totalRho[j][k];
	      res += sqr(a1_local[j][k]*phi[j-1][k]
			 +b1_orig[j][k]*phi[j][k]
			 +c1_local[j][k]*phi[j+1][k]+
			 a2_local[j][k]*phi[j][k-1]
			 +b2_orig[j][k]*phi[j][k]+
			 totalRho[j][k]);

	      foo[j][k] = sqr(a1_local[j][k]*phi[j-1][k]
			      +b1_orig[j][k]*phi[j][k]
			      +c1_local[j][k]*phi[j+1][k]+
			      a2_local[j][k]*phi[j][k-1]
			      +b2_orig[j][k]*phi[j][k]+
			      totalRho[j][k]);
	    }
	}

      for(k=1;k<K;k++)
	{
	  j = 0;
	  if (b1_orig[j][k] !=0)
	    {
	      temp += totalRho[j][k];
	      res += sqr(b1_orig[j][k]*phi[j][k]+c1_local[j][k]*phi[j+1][k]+
			 a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]+
			 c2_local[j][k]*phi[j][k+1]+
			 totalRho[j][k]);

	      foo[j][k] = sqr(b1_orig[j][k]*phi[j][k]
			      +c1_local[j][k]*phi[j+1][k]+
			      a2_local[j][k]*phi[j][k-1]
			      +b2_orig[j][k]*phi[j][k]
			      +c2_local[j][k]*phi[j][k+1]+
			      totalRho[j][k]);
	    }
	  
	  j = J;
	  if (b1_orig[j][k] !=0)
	    {
	      temp += totalRho[j][k];
	      res += sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]+
			 a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]+
			 c2_local[j][k]*phi[j][k+1]+
			 totalRho[j][k]);
				
	      foo[j][k] = sqr(a1_local[j][k]*phi[j-1][k]
			      +b1_orig[j][k]*phi[j][k]+
			      a2_local[j][k]*phi[j][k-1]
			      +b2_orig[j][k]*phi[j][k]
			      +c2_local[j][k]*phi[j][k+1]+
			      totalRho[j][k]);
	    }
	}
      
      // interior
      for (j=1;j<J;j++) for(k=1;k<K;k++)
	{		
	  res += sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]
		     +c1_local[j][k]*phi[j+1][k]+
		     a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]
		     +c2_local[j][k]*phi[j][k+1]+
		     totalRho[j][k]);

	  foo[j][k] = sqr(a1_local[j][k]*phi[j-1][k]+b1_orig[j][k]*phi[j][k]
			  +c1_local[j][k]*phi[j+1][k]+
			  a2_local[j][k]*phi[j][k-1]+b2_orig[j][k]*phi[j][k]
			  +c2_local[j][k]*phi[j][k+1]+
			  totalRho[j][k]);
			
	  
	  if (b1_orig[j][k]!=0)
	    temp += sqr(totalRho[j][k]);
	}
      if (temp)
	{
	  res = sqrt(res/temp);
	  for (j=0;j<=J;j++) for(k=0;k<=K;k++)
	    foo[j][k] = sqrt(foo[j][k]/temp);
	  
	}
		
      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:
      for (j=0;j<=J;j++)
	for (k=0;k<=K;k++)
	  resCoeff_local[j][k][4] = resCoeff[j][k][4];
		
      res = psolve->Resid(totalRho, phi);	
    }	
  return (res);
}

Vector2 Boltzmann::ParticleFlux(Scalar dt)
{
  Scalar PartCurrent;
  Scalar EnergyCurrent;
  Scalar current = 0.;
  int j,k;

  Boundary* bPtr;
  Boundary*** bc1 = grid->GetBC1();
  Boundary*** bc2 = grid->GetBC2();
  Vector3** dS = grid->dS();

  const Scalar scale = dt/BoltzSpecies->get_q()*.5*sqrt(2*temperature*fabs(ELECTRON_CHARGE)/M_PI/BoltzSpecies->get_m());

  PartCurrent = 0;
  EnergyCurrent = 0;