pdc1.c

一维等离子体静电粒子模拟程序

/**********************************************************************
**	PDC1 - ELECTROSTATIC 1 DIMENSIONAL BOUNDED PLASMA SIMULATION
**      Simulates an electrostatic plasma in cylindrical geometry.
**	Includes an external LCR circuit, and optionally a magnetic field along
**      the z-axis. 
**	
**	REVISION/PROGRAMMER/DATE
**	
**	1.0 modified version of PDP1 in cylindrical geometry
**      /M. V. Alves and Vahid Vahedi/07-15-89
**
**	2.0 Addition of dumping and restoring a state of the system
**	/Vahid Vahedi and M. V. Alves/05-31-90
**
**      2.1 Conversion from IBM PC to UNIX 
**      /Payam Mirrashidi/
**     
**      2.2 11/28/94 New Xgraphics; Couple of minor bug fixes. P.G.
**      Lotsa new disgnostics, attempt to make it like pdp1
**
**      3.0 2/22/99 Programs extensively re-written and re-organised
**         /Helen Smith/
**         New features include:  
**      1) calculations all use normalised parameters and only    
**         diagnostics are re-normalised  (to reduce calc time)   
**      2) A ramp is now available for the ac voltage             
**      3) Simulation particle numbers can be re-scaled to keep   
**         them within specified parameters                       
**	4) New diagnostics, such as sheath width, plasma impedance and
**	   angle btwn V-I 
**      5) Includes more gases.
**
**      3.1 March 2000
**      Code modified to allow modelling of cylindrical glow discharge
**      with no central conductor.
**      1) particle mover changed. Uses rotation in cylindrical coords
**      (see B&L pg 338) 
**      2) a new poisson solution was added for r0 = 0 
**      3) the end cell volumes were adjusted, in order to give the correct
**         weightings.
**      4) substepping for the ions was included 
**
**      3.1 March 2000
**      Code modified to allow modelling of cylindrical glow discharge
**      with no central conductor.
**      1) particle mover changed. Uses rotation in cylindrical coords
**      (see B&L pg 338)
**      2) a new poisson solution was added for r0 = 0
**      3) the end cell volumes were adjusted, in order to give the correct
**         weightings.
**      4) substepping for the ions was included
**
**     June 2000
**     1) Some small bug fixes were made: yp(0) now uses adjusted cell volume,
**        np_start is made an integer value.
**     2) Non-uniform grid spacings included - code will now run uniform or
**        non-uniform depending on value set for parameter unif set in input
**        file  (1=uniform 0 = non-uniform). HOWEVER there is still a problem
**        with the non-uniform grid which is probably because the Poisson eqn
**        is only to order dt. Need to derive 2nd order Poisson soln for
**        better accuracy.
**     3) Particle densities can now be loaded with several profiles,
**        by specifying the value of "profile" read in from the species
**        data - 0) uniform, 1) cosine, 2) J0. # particles to be loaded
**        now depends on density, reactor volume and profile.
**     4) velocity diagnostic commented out, since it takes ~ 40% of cpu
**        time and is not commonly considered. Might need to add a flag
**        if it is required later.
**
**     Still problems with rescale routine, so re-scaling currently turned off. 
**     
**     June 2000
**       /Hae June Lee/
**     Radiation transport model has been included. The RT part is mainly
**     implanted in the radtr.c.
**
**     June 2000 Addition of Parallelization scheme.
**       /Emi Kawamura/
**       The parallelization scheme is documented in UCB/ERL Report No. M99/58.
**       The parallel code requires the use of the MPI libraries available at 
**       http://www.mcs.anl.gov/mpi. In order to turn on the parallel features
**       the executable must be compiled with the -DMPI_1D flag. Please
**       also read the file README.parallel.

***********************************************************************/ 
#include "pdc1.h"
#include <xgrafix.h>

void display_title(void);
void InitWindows(void);
void XGMainLoop(void);
void rescale(int);

#ifdef POW
#define POW_DT 10000
void pow_balance(void);
char powfile[50];
float *pow_bal, *powbal_tot, pow_sum;
int n_pow;
#endif

int main(int argc, char *argv[])
{
  int i, isp;

#ifdef MPI_1D
  MPI_Init(&argc, &argv);                    /* Initiate MPI            */
  MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);   /* Get my process rank     */
  MPI_Comm_size(MPI_COMM_WORLD, &num_procs); /* Get number of processes */
  if (my_rank == ROOT) printf("Using Parallel Version of PDC1\n");
  startwtime = MPI_Wtime();
#else
  startwtime = clock();
#endif

  /* Set min number of particles for simulation */
  min_np = 500;     /* minimum # particles before rescaling */
 
#ifdef MPI_1D
  if (my_rank == ROOT) display_title();      /*	Display XPDC1 title 	 */
#else
  display_title();                           /* Display XPDC1 title      */
#endif

  XGInit(argc,argv,&atime); /* Initialize XGrafix stuff      */

#ifdef MPI_1D
  if(!WasDumpFileGiven) {
    if (my_rank == ROOT) printf("Parallel version needs dumpfile,
      please generate one using non-parallel XPDC1\n");
    MPI_Finalize();
    exit(1);
  }
  if(theRunWithXFlag) {
    if (my_rank == ROOT) printf("Parallel version should be run
      without XGrafix. Please use -nox flag.\n");
    MPI_Finalize();
    exit(1);
  }
#endif

  start();   
  if(theRunWithXFlag) 
  {
     mccdiag_init();  /* Initialize MCC rate diagnostic*/
     
   #ifdef MPI_1D
     if (my_rank == ROOT) InitWindows();  /* Initialize diagnostic windows   */
   #else
     InitWindows();  /* Initialize diagnostic windows   */
   #endif
  }

  init_rho();

#ifdef MPI_1D
  /* Sum densities of procs to get total density and send result to all CPUs */
  for (isp = 0; isp < nsp; isp++) {
    MPI_Allreduce(sp_n[isp], sp_n_tot[isp], ng, MPI_FLOAT, MPI_SUM, 
      MPI_COMM_WORLD);
      
    for (i=0;i<ng;i++) sp_n_tot[isp][i] /= vol[i];
  }

  /* Print out total density of first two particle species.*/
  if (my_rank == ROOT) {
    printf("sp_n_tot[0][ng/2] = %g\n", sp_n_tot[0][ng/2]);
    printf("sp_n_tot[1][ng/2] = %g\n", sp_n_tot[1][ng/2]);
    printf("proc 0: sp_n[0][ng/2] = %g\n", sp_n[0][ng/2]/vol[ng/2]);
    printf("proc 0: sp_n[1][ng/2] = %g\n", sp_n[1][ng/2]/vol[ng/2]);
  }
#else
  for (isp = 0; isp < nsp; isp++) 
    for (i=0;i<ng;i++) sp_n[isp][i] /= vol[i];
   printf("sp_n[0][ng/2] = %g\n", sp_n[0][ng/2]);
   printf("sp_n[1][ng/2] = %g\n", sp_n[1][ng/2]);
#endif

  field_init();
  fields();            	/* Initialize potential & fields  */
  
  if(theRunWithXFlag) 
      history();    /* Initialize time histories     */

#ifdef MPI_1D
  time0 = MPI_Wtime();
#else
  time0 = clock();
#endif

#ifdef POW
  strcpy(powfile, theDumpFile);
  strcat(powfile, ".pow");
  n_pow = 2*nsp + 4;
  pow_bal = (float *)malloc(n_pow*sizeof(float));
  powbal_tot = (float *)malloc(n_pow*sizeof(float));
  pow_sum = 0.0;
  for (i=0; i<n_pow; i++) pow_bal[i] = powbal_tot[i] = 0.0;
#endif

  XGStart();		/* Start XGrafix main loop       */
  
#ifdef MPI_1D
  if (my_rank == ROOT)
    for(i=0; i<nsp; i++) printf("proc 0: np[%d] = %d, \n", i, np[i]);
#else
  for(i=0; i<nsp; i++) printf("np[%d] = %d, \n", i, np[i]);
#endif

  return (0);
}/* End of MAIN */


/********************************************************/
/*	The main physics loop				*/

void XGMainLoop(void)
{   
    int isp, j, k;
    float frac;

    atime += dt; it_rate++;

/* loop over particle species */
    for (isp=0; isp<nsp; isp++)
    {
/* if # sub-steps == main step then move particles */
       if (!(k_count[isp]%sp_k[isp]))
       {
         it[isp]++;             /* Advance by one time-step */
         move(isp);             /* Advance particle positions and velocities */
         wall(isp);             /* Remove particles that hit the walls */
         (*mccptr) (isp);       /* Carry out MCC collisions */
         gather(isp);           /* Calculate new potentials and fields */
         k_count[isp]=0;
      }
      k_count[isp]++;
      frac = ((float)k_count[isp])/sp_k[isp];
      for (j=0; j<ng; j++)
         sp_n[isp][j] = (1 - frac)*sp_n_0[isp][j] + frac*sp_n_k[isp][j] 
                       + sp_n_mcc[isp][j];

    #ifdef MPI_1D
      /* Sum N(x) of CPUs to get tot. density and send result to all CPUs*/
      MPI_Allreduce(sp_n[isp], sp_n_tot[isp], ng, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
      for (j=0;j<ng;j++) sp_n_tot[isp][j] /= vol[j];
    #else
      for (j=0;j<ng;j++) sp_n[isp][j] /= vol[j];
    #endif

    }

  #ifdef MPI_1D
    /* Sum iwall of CPUs to get tot. iwall and send result to all CPUs*/
    MPI_Allreduce(iwall, iwall_tot, NSMAX, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
  #endif 
    
    fields();             /* Calculate new potentials and fields */

    /*-------------- Evoultion of excited states ----------------*/
    if (RT_flag) {
      if (ex_count%sp_k_ex==0) {
        for (j=0;j<nc_RT;j++) eden_ave[j]=0.5*eden_ave[j]/sp_k_ex;
        excited_evolution();                     /* LEEHJ */
        ex_count=0;
      }
      ex_count++;
      
      if (unif) for (j=0;j<nc_RT;j++) {
        isp=j*nc_ratio;
        for (k=0;k<nc_ratio;k++) {
	#ifdef MPI_1D
           eden_ave[j]+=sp_n_tot[0][isp]+sp_n_tot[0][isp+1];
        #else
           eden_ave[j]+=sp_n[0][isp]+sp_n[0][isp+1];
        #endif
           isp++;
        }
      }
    }
    if(theRunWithXFlag) history();

 #ifdef POW
    pow_balance();
 #endif

  /* Check to see if particle numbers are within bounds, and rescale if not */
  /* Rescaling currently turned off - needs to be debugged */
    
  for (isp=0; isp<nsp; isp++)
    /*
    if ((np[isp] > 0.9*maxnp[isp]) || (np[isp] < min_np)) rescale(isp);
    
    if ((np[isp] > 0.9*maxnp[isp]) || (np[isp] < min_np))
    {
       
  #ifdef MPI_1D
     if (my_rank == ROOT) 
     {
     if(np[isp] > 0.9*maxnp[isp]) { printf("Too many particles - stopping simulation \n"); }
     else  { printf("Too few particles - stopping simulation \n"); }
     }
  #else
  if(np[isp] > 0.9*maxnp[isp]) { printf("Too many particles - stopping simulation \n"); }
     else  { printf("Too few particles - stopping simulation \n"); }
  #endif
     getchar();
     exit(1);
    }
    */

/* Print out diagnostics every 5000 electron time-steps */
#ifdef MPI_1D
   if (my_rank == ROOT && it[0]%5000 == 0)
     printf("For proc 0: t = %e it = %d n_e = %d n_i = %d Ez = %f Iz = %f\n",atime,it[0],np[0],np[1], E_z/(ptopn*dr), I_z);
#else
   if (it[0]%5000 == 0)
     printf("t = %e it = %d n_e = %d n_i = %d Ez = %f Iz = %f\n",atime,it[0],np[0],np[1], E_z/(ptopn*dr), I_z);
#endif
}	
/**************************************************************/
void display_title()
{
  printf("\n\nPDC1 - Cylindrical Bounded Electrostatic 1 Dimensional Code\n");
  printf("version 3.1\n");
  printf("(c) Copyright 1989-95 Regents of the University of California\n");
  printf("Plasma Theory and Simulation Group\n");
  printf("University of California - Berkeley\n");
  puts("See README file for information about Licensing\n");
}

/*****************************************************************/
#ifdef POW
void pow_balance()
{
 int i;
 FILE *PowFile;
 float coll_loss, wall_loss, temp;
 
 if (it[0]%POW_DT == 0) {
 #ifdef MPI_1D
   
   MPI_Allreduce(pow_bal, powbal_tot, n_pow, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
   if (my_rank == ROOT) {

     PowFile = fopen(powfile, "a");
     fprintf(PowFile, "time = %e\n", atime);
     for(i=0; i<nsp; i++) {
       fprintf(PowFile, "species=%d, particle loss to wall =%g\n", 
               i, nc2p*powbal_tot[i]/sp_k[i]);
     }
     fprintf(PowFile, "number of  ionization events =%g\n",
             nc2p*powbal_tot[2*nsp]);
     
     wall_loss = 0.0;
     coll_loss = 0.0;
     for(i=0; i<nsp; i++) {
       temp = nc2p*ECHARGE*powbal_tot[nsp+i]/(POW_DT*dt*sp_k[i]);
       wall_loss += temp;
       fprintf(PowFile, "species=%d, power loss to wall = %g W\n", i, temp);
     }
   
     temp = nc2p*ECHARGE*(eion*powbal_tot[2*nsp] + eexc*powbal_tot[2*nsp+1])/
            (POW_DT*dt);
     coll_loss += temp;
     fprintf(PowFile, "power loss to inelastic colls = %g W\n", temp);
     
     temp = nc2p*ECHARGE*powbal_tot[2*nsp+2]/(POW_DT*dt);  
     coll_loss += temp;
     fprintf(PowFile, "power loss to elastic colls = %g W\n", temp);
   
     temp = nc2p*ECHARGE*powbal_tot[2*nsp+3]/(POW_DT*dt*sp_k[1]);
     coll_loss += temp;
     fprintf(PowFile, "power loss to charge exchange colls = %g W\n", temp);
   
     fprintf(PowFile, "total power loss to colls = %g W\n", coll_loss);
     fprintf(PowFile, "total power loss to wall = %g W\n", wall_loss);
     fprintf(PowFile, "total power loss = %g W\n", coll_loss+wall_loss);
     fprintf(PowFile, "power Input = %g W\n", pow_sum/POW_DT);
     fprintf(PowFile, "\n");
     fclose(PowFile);
   }
 #else
   PowFile = fopen(powfile, "a");
   fprintf(PowFile, "time = %e\n", atime);
   for(i=0; i<nsp; i++) {
     fprintf(PowFile, "species=%d, particle loss to wall =%g\n",
             i, nc2p*pow_bal[i]/sp_k[i]);
   }
   fprintf(PowFile, "number of ionization events =%g\n", nc2p*pow_bal[2*nsp]);
     
   wall_loss = 0.0;
   coll_loss = 0.0;
   for(i=0; i<nsp; i++) {
     temp = nc2p*ECHARGE*pow_bal[nsp+i]/(POW_DT*dt*sp_k[i]);
     wall_loss += temp;
     fprintf(PowFile, "species=%d, power loss to wall = %g W\n", i, temp);
   }
   
   temp = nc2p*ECHARGE*(eion*pow_bal[2*nsp]+eexc*pow_bal[2*nsp+1])/(POW_DT*dt);
   coll_loss += temp;
   fprintf(PowFile, "power loss to inelastic colls = %g W\n", temp);
   
   temp = nc2p*ECHARGE*pow_bal[2*nsp+2]/(POW_DT*dt);  
   coll_loss += temp;
   fprintf(PowFile, "power loss to elastic colls = %g W\n", temp);
   
   temp = nc2p*ECHARGE*pow_bal[2*nsp+3]/(POW_DT*dt*sp_k[1]);
   coll_loss += temp;
   fprintf(PowFile, "power loss to charge exchange colls = %g W\n", temp);
   
   fprintf(PowFile, "total power loss to colls = %g W\n", coll_loss);
   fprintf(PowFile, "total power loss to wall = %g W\n", wall_loss);
   fprintf(PowFile, "total power loss = %g W\n", coll_loss+wall_loss);
   fprintf(PowFile, "power Input = %f W\n", pow_sum/POW_DT);
   fprintf(PowFile, "\n");
   fclose(PowFile);
 #endif
   for(i=0; i<n_pow; i++) pow_bal[i] = 0.0;
   pow_sum = 0.0;
 }
 else {
   pow_sum += E_z*I_z*height/(ptopn*dr);
   
   for (i=0; i<nsp; i++) pow_bal[i] += n_lost[i];
   for (i=0; i<nsp; i++) pow_bal[nsp+i] += e_lost[i];
   pow_bal[2*nsp]   += nu_iz;
   pow_bal[2*nsp+1] += nu_ex;
   pow_bal[2*nsp+2] += els_loss;
   pow_bal[2*nsp+3] += cx_loss;
 }
   
} 
#endif