start.c

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

/**********************************************************************
**   start
**   Read in input file, allocate arrays, normalise variables,
**   restore data from a previous run or load in a new particle
**   distribution and various other "setting up" steps
**   Parameter grid determines type of grid to be loaded:
**   grid = 0 non-uniform (constant volume) mesh
**   grid = 1 uniform mesh
**   grid = 2 linearly decreasing mesh
**   grid = 3 uniform mesh with 2 grid spacings
***********************************************************************/
#include "pdc1.h"
#include <xgrafix.h>

FILE *InputDeck;

void species(int initnp[NSMAX][2], int isp);
void Restore(char *filename);
int grid, nc0, profile;
float iinitn, dr1;
void unif_load(int initnp[NSMAX][2]);
void cos_load(int initnp[NSMAX][2], float iinitn);
void bes_load(int initnp[NSMAX][2], float iinitn);
void input_for_rad();

/*------------------------------------------------------------*/
void start(void)
{
  float Bz, f0;
  float ntemp, sum_nc, grad;
  char aachar[100];
  int i, j, isp, nc1, initnp[NSMAX][2];
  
/*---------------------------------------------*/
/* Open files and read in "general" parameters */
  
  if (!WasInputFileGiven) 
  {
  #ifdef MPI_1D
    if(my_rank == ROOT){
      puts("Syntax:");
      puts("\tPDC1 -i <inputdeck> -d [dumpfile]");
    }
  #else
    puts("Syntax:");
    puts("\tPDC1 -i <inputdeck> -d [dumpfile]");
  #endif
    exit(1);
  }
  InputDeck = fopen(theInputFile, "r");

  /*---------------------------------------------*/
  /* Read lines until we get to numbers */
  while (fscanf(InputDeck,"%d %d %d %f %f %f %f %f %f %f",
		&nsp, &nc, &grid, &nc2p, &dt, &r0, &r1, &height, &epsilon, &Bz) <10)
    fscanf(InputDeck, "%s", aachar);
  
  while (fscanf(InputDeck, "%f %f %f %f %f %f %f",
		&rhoback, &backgrnd_i, &dde, &extr, &extl, &extc, &q_0) < 7)
    fscanf(InputDeck, "%s", aachar);
  
  while (fscanf(InputDeck, "%d %c %f %f %f %f %f %f", &ramped,
		&src, &dcbias, &dc_ramp, &acbias, &ac_ramp, &f0, &theta0) <8) 
    fscanf(InputDeck, "%s", aachar);

 if((src == 'W') ||(src == 'K'))
    {
      while (fscanf(InputDeck, "%d %f %f %f %f %f %f", &ramped2,
                &dcbias2, &dc_ramp2, &acbias2, &ac_ramp2, &f02, &theta02) <7) 
        fscanf(InputDeck, "%s", aachar);
    } 

  
  while (fscanf(InputDeck, "%d %d %d %d %d %d %d", &secondary, &ecollisional,
		&icollisional, &reflux, &nfft, &nsmoothing, &RT_flag) <7)
    fscanf(InputDeck, "%s", aachar);
  
  while (fscanf(InputDeck, "%f %f %d %f %f %d", &seec[0], &seec[1],
		&ionspecies, &pressure, &T_gas, &gas) <6)
    fscanf(InputDeck, "%s", aachar);
  
  /*---------------------------------------------*/
  /* Check for errors in the "general" input parameters */
  
  #ifdef MPI_1D
  if (RT_flag) {
    printf("For parallel version, RT simulation is not possible\n");
    exit(1);
  }
  #endif
  if (extl < 0. || extr < 0. || extc < 0.)
  {
  #ifdef MPI_1D
    if(my_rank == ROOT)
      printf("START-ERROR- external circuit improperly defined R,L,C < 0\n");
  #else
    printf("START-ERROR- external circuit improperly defined R,L,C < 0\n");
  #endif
    exit(1);
  }
  if (!dt)
  {
  #ifdef MPI_1D
    if(my_rank == ROOT)
      printf("START: dt == 0!!\n");
  #else
    printf("START: dt == 0!!\n");
  #endif
    exit(1);
  }
  if (nsp > NSMAX)
  {
  #ifdef MPI_1D
    if(my_rank == ROOT)
      printf("START: nsp > NSMAX.\n");
  #else
    printf("START: nsp > NSMAX.\n");
  #endif
    exit(1);
  }
    
  /* Check to see if nfft is an integer power of 2 */
  if(nfft)
  {
    for(i=0; i<= 31; i++)
    {
      if(nfft == (1<<i)) 
      {
	i=0;
	break;
      }
    }
    if(i)
    {
    #ifdef MPI_1D
      if(my_rank == ROOT)
        printf("START: nfft is not an integer power of 2.\n");
    #else
      printf("START: nfft is not an integer power of 2.\n");
    #endif
      exit(1);
    }
  }

  ng= nc+1;
  
  /*-----------------------------------*/
  /*  Set collision parameters         */
   mccptr =argonmcc;         /* gas uses argonmcc (without action) if collisions are off */
   if (ecollisional || icollisional)
   {
      switch (gas)
      {
         case ARGON:
	   mccptr =argonmcc; 
	   sigmaMptr =asigma1a; 
	   sigmaTptr =asigmaT; 
	   mueFacPtr = mue_factor; 
	   eion = 15.76;
	   eexc = 11.62;
	   e_me = 11.55;
	   e_miz = 4.16;
        #ifdef MPI_1D
	   if(my_rank == ROOT)
	     printf("Using anisotropic HBS cross-section set\n");
	#else
           printf("Using anisotropic HBS cross-section set\n");
        #endif
	   break;

         case ARGON_LK:
	   mccptr =argonmccLK; 
	   sigmaMptr =LK_asigma1a; 
	   sigmaTptr =LK_asigmaT; 
	   mueFacPtr = mue_factorLK; 
	   eion = 15.76;
	   eexc = 11.62;
	   e_me = 11.55;
	   e_miz = 4.16;
        #ifdef MPI_1D
	   if(my_rank == ROOT)
	     printf("Using isotropic L&K cross-section set\n");
	#else
           printf("Using isotropic L&K cross-section set\n");
        #endif
	   break;
      }
   }  


  /*---------------------------------------------*/
  /* Calculate normalised cell radii and volumes */

  r_grid = (float  *)malloc(ng*sizeof(float));
  dr_n = (float  *)malloc(ng*sizeof(float));

/* Constant cell volume mesh */
  if (grid == 0)
  {
      j_pointer = const_vol_grid;

      dr = -r0 + sqrt(r0*r0 + (r1*r1 - r0*r0)/nc);
      r0 /= dr;
      r1 /= dr;

      r_grid[0] = r0;
      r_grid[1] = r0+1;
      for (i=2; i<= nc; i++) r_grid[i] = sqrt(2.0*r_grid[i-1]*r_grid[i-1] - r_grid[i-2]*r_grid[i-2]);

      for (i=0; i< nc; i++) dr_n[i] = r_grid[i+1] - r_grid[i];
  }
/* Constant uniform grid */
  else if(grid == 1)
  {
      j_pointer = unif_grid;
      dr = (r1-r0)/nc;
      r0 /=dr;
      r1 /=dr;

      for ( i=0; i<= nc; i++) r_grid[i] = r0 + i;
      for ( i=0; i<= nc; i++) dr_n[i] = 1;
  }
/* Linear dr mesh */
  else if(grid == 2)
  {
      j_pointer = linear_grid;

/* determine first dr */
      dr = 1.5*(r1 - r0)/nc; /* choose the first grid cell spacing - between 1.2 and 2.0 works best */
      r0 /= dr;
      r1 /= dr;

      sum_nc = 0.0;
      for (i=1; i<nc; i++) sum_nc += i;
      grad = (r1 - r0 - nc)/sum_nc;
      for (i=0; i< nc; i++) dr_n[i] = grad*i + 1.0;
      r_grid[0] = r0;
      for (i=1; i<= nc; i++) r_grid[i] = r_grid[i-1] + dr_n[i-1];
  }
/* Uniform mesh with two spacings */
  else if(grid == 3)
  {
     j_pointer = split_grid;

     nc0 = 0.2*nc;
     nc1 = nc - nc0;
     dr = 0.5*(r1 - r0)/nc0;                      /* 1/4 length with nc0 cells */
     dr1 = 0.5*(r1 - r0)/nc1;                     /* 3/4 length with nc1 cells */
     r0 /= dr;
     r1 /= dr;

     for (i=0; i< nc0; i++) dr_n[i] = 1.0;         /* normalise by dr */
     for (i=nc0; i< nc; i++) dr_n[i] = dr1/dr;
     r_grid[0] = r0;
     for (i=1; i<= nc; i++) r_grid[i] = r_grid[i-1] + dr_n[i-1];
  }

  vol = (float  *)malloc(ng*sizeof(float));
  vol[0] = r0 + 1.0/3.0; 
  vol[nc] = dr_n[nc-1]*(r_grid[nc] - dr_n[nc-1]/3.0); 
  for ( i=1; i< nc; i++) 
    vol[i] = r_grid[i]*(dr_n[i] + dr_n[i-1]) + (dr_n[i]*dr_n[i] - dr_n[i-1]*dr_n[i-1])/3.0;

  /*------------------------------------*/
  /* Conversion to dimensionless units */
  
  gas_den = NperTORR*pressure/(T_gas +DBL_MIN);
  area0 = TWOPI*r0*dr*height;
  theta0 *= TWOPI/360;
  Trf =1/(f0+DBL_MIN);
  cycle = Trf/dt;
  w0 = TWOPI*f0;
  w02 = TWOPI*f02;
  epsilon*= EPS0;
  q_1 = q_2 = q_3 = q_0;
  reactor_vol = PI*height*dr*dr*(r1*r1 - r0*r0);  /* in real units */
  rescale_on=1;
  
  /*-----------------------------------------------*/
  /* Calculate normalisation parameters */
  me_e = EMASS/ECHARGE;
  vscale = dr/dt;
  ptopn = (ECHARGE/EMASS)*dt*dt/(dr*dr);  /* Normalisation for potential */
  dntod =  nc2p/(PI*height*dr*dr); /* Normalisation for density */

  /* Check ramped source parameters                      */
  if(ramped)
  {
    if(dc_ramp*dcbias < 0.0)
    {
     #ifdef MPI_1D
       if(my_rank == ROOT){
         fprintf(stderr, "\nWARNING: Ramp and dcbias have opposite signs.");
         fprintf(stderr, "\nThe external source will be discontineous!\n\n");
       }
     #else
       fprintf(stderr, "\nWARNING: Ramp and dcbias have opposite signs.");
       fprintf(stderr, "\nThe external source will be discontineous!\n\n");
     #endif
    }
    if (dc_ramp) risetime= fabs(dcbias)/fabs(dc_ramp);
  }
 if(ramped2)
  {
    if(dc_ramp2*dcbias2 < 0.0)
    {
     #ifdef MPI_1D
       if(my_rank == ROOT){
         fprintf(stderr, "\nWARNING: Ramp and dcbias have opposite signs.");
         fprintf(stderr, "\nThe external source will be discontineous!\n\n");
       }
     #else
       fprintf(stderr, "\nWARNING: Ramp and dcbias have opposite signs.");
       fprintf(stderr, "\nThe external source will be discontineous!\n\n");
     #endif
    }
    if (dc_ramp2) risetime2= fabs(dcbias2)/fabs(dc_ramp2);    
  }


  /*-----------------------------------*/
  /* Allocate space for field parameters */
  
  r_array= (float  *)malloc(ng*sizeof(float));
  for(i=0; i< ng; i++) r_array[i]= r_grid[i]*dr;

  sp_n= (float  **)malloc(nsp*sizeof(float  *));
  sp_n_k= (float  **)malloc(nsp*sizeof(float  *));
  sp_n_0= (float  **)malloc(nsp*sizeof(float  *));
  sp_n_mcc= (float  **)malloc(nsp*sizeof(float  *));
  for(i=0; i<nsp; i++)