argonmcc_hbs.c

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

/* Monte-Carlo collision routines for particle species with background gas and
** each other.
** asigma1:  Ar + e -> Ar + e      scattering
** asigma2:  Ar + e -> Ar* + e     lumped excitation cross-section (inc 2 radiative and 2 metastable levels)
** asigma3:  Ar + e -> Ar+ + 2e     ionisation from ground state
** asigma4:  Ar+ + Ar -> Ar + Ar+   charge transfer (ions)
** asigma5:  Ar+ + Ar -> Ar+ + Ar   elastic scattering of ions
** asigma6:  Ar^m + e -> Ar+ + 2e   Ionization of excited states by electron impact
** asigma7:  Ar + e -> Ar^m + 2e   Excitation to metastable states by electron impact (subset of asigma2)
*/

#include "pdc1.h"
#include "xgrafix.h"

#define NEMAX    100
#define NIMAX    100

void newvel(float, float, float *, float *, float *, int);

/* Set factor for calculating electron mobility */
float mue_factor()
{
  return(1.25); /* anisotropic HBS cs set */
}

/*-------------------------------------------------------------*/
/* Monte Carlo scheme for electron and ion-neutrals collisions */

void argonmcc(int isp)
{
  static int init_flag=1, icorrect=0;
  static float ecol_extra, icol_extra, ecol_extra2;		
  static float max_sigmav_e, max_sigmav_i, max_sigmav_e2;	
  static float col_prob_e, col_prob_i, Cdt, N_ex_max=1;
  
  register int j, i;
  int N, k, nnp, s, index;
  float etov, dt_correct;
  float rand_num, v_sqrd, vel, sigma_total, sum_prob;
  float engy, rengy, frac, phi1;
  float cosphi, sinphi, coschi, sinchi;
  float vneutr, vneutth, vneutz;
  float temp;						

  if (init_flag)
  {
    /*---------------------------------------*/
    /* Calculating the null collision prob.  */
    
    if (ionspecies) 
    {
       ionsp = 1;
       vg_therm = sqrt(T_gas/(Escale[ionsp]));
    }
    if(ecollisional)
    {
      ecolsp = 0;
      max_sigmav_e = 0.0;
      etov = 2.0*ECHARGE/mass[ecolsp];
      for (i=0; i< 10*NEMAX; i++)
      {
	engy = 0.1*i;
	max_sigmav_e= max(max_sigmav_e, sqrt(engy*etov)*(asigma1(engy)+
                      asigma2(engy)+asigma3(engy)));
      }
      col_prob_e = 1.0 -exp(-max_sigmav_e*gas_den*sp_k[ecolsp]*dt);

  #ifdef MPI_1D
      if (my_rank == ROOT){
        printf("Gas density = %g\n", gas_den);
        printf("Electron coll prob = %g, max(sigma*v) = %g, nsigmavdt= %g\n", 
            col_prob_e, max_sigmav_e, max_sigmav_e*gas_den*dt*sp_k[ecolsp]);
      }
  #else                      
      printf("Gas density = %g\n", gas_den);
      printf("Electron coll prob = %g, max(sigma*v) = %g, nsigmavdt= %g\n", 
            col_prob_e, max_sigmav_e, max_sigmav_e*gas_den*dt*sp_k[ecolsp]);
  #endif

      max_sigmav_e2 = 0.0;		
/* If radiation transport turned on calculate ionisation of metastables by
electron impact */
      if (RT_flag) 
      {
        for (i=0; i<= 2000; i++) 
        {
          engy = 0.1*i;
          max_sigmav_e2= max(max_sigmav_e2, sqrt(engy*etov)*asigma6(engy));
        }
        Cdt=max_sigmav_e2*dt*sp_k[ecolsp];
      }
    }
    
    if(icollisional)
    {
      icolsp = 1;
      max_sigmav_i = 0.0;
      etov = 2.0*ECHARGE/mass[icolsp];
      for (i=0; i<10*NIMAX; i++)
      {
	engy = 0.1*i;
	max_sigmav_i= max(max_sigmav_i, sqrt(engy*etov)
			  *(asigma4(engy)+asigma5(engy)));
      }
      col_prob_i = 1 -exp(-max_sigmav_i*gas_den*sp_k[icolsp]*dt);

  #ifdef MPI_1D
      if (my_rank == ROOT) /* MPI modification; only ROOT prints */
        printf("Ar+ coll prob = %g, max(sigma*v) = %g, nsigmavdt= %g\n", col_prob_i, max_sigmav_i, max_sigmav_i*gas_den*dt*sp_k[icolsp]);
  #else
      printf("Ar+ coll prob = %g, max(sigma*v) = %g, nsigmavdt= %g\n", col_prob_i, max_sigmav_i, max_sigmav_i*gas_den*dt*sp_k[icolsp]);
  #endif

    }
    ecol_extra = icol_extra = ecol_extra2 = 0.5;
    init_flag = 0;
  }
  
  /*------------------------------------------------------------* 
       Electron Collision : Ar* + e  -> Ar+ + 2e  ( by LEEHJ )
   *------------------------------------------------------------*/
  if(RT_flag && ecollisional && isp==ecolsp) 
  {
    dt_correct=gas_den/(N_ex_max*10);
    icorrect++;
    for (N_ex_max=0,j=0;j<nc_RT;j++) N_ex_max=max(N_ex_max, sp_ex[j]+sp_exm[j]);
    if (dt_correct<=icorrect) 
    {
      if(theRunWithXFlag) for (j=0; j<ng; j++) mccrate[5][j] = 0;
      ecol_extra2 += np[ecolsp]* (1-exp(-N_ex_max*icorrect*Cdt));
      N=ecol_extra2;
      ecol_extra2-=N;
  
      nnp = np[ecolsp];
      for(j=0; j< N; j++)
      {
        nnp--;
        index= nnp*frand();
  
        swap(r[ecolsp],index,nnp);
        swap(v_r[ecolsp],index,nnp);
        swap(v_theta[ecolsp],index,nnp);
        swap(v_z[ecolsp],index,nnp);
  
        v_sqrd = v_r[ecolsp][nnp]*v_r[ecolsp][nnp] 
               + v_theta[ecolsp][nnp]*v_theta[ecolsp][nnp]  
               + v_z[ecolsp][nnp]*v_z[ecolsp][nnp];
        engy= Escale[ecolsp]*v_sqrd;
        vel = sqrt(v_sqrd);
        sigma_total = max_sigmav_e2/(vel*vscale);
        rand_num= frand();     

        s=j_pointer(ecolsp,nnp)/nc_ratio; 
        temp=(sp_ex[s]+sp_exm[s])/N_ex_max;
        if (engy>e_miz && rand_num <= (sum_prob=asigma6(engy)/sigma_total)*temp)
        {
        /* first normalize the velocity components by the speed, vel */
          v_r[ecolsp][nnp] /= vel;
          v_theta[ecolsp][nnp] /= vel;
          v_z[ecolsp][nnp] /= vel;

          engy -= e_miz;
          rengy = 2.0*tan(frand()*atan(engy/20.0));
          engy -= rengy;
      
          /*-----------------------------------*/
          /* scatter the created electron	    */

          vel = sqrt(rengy/Escale[ecolsp]);
          k = np[ecolsp];
          v_r[ecolsp][k] = v_r[ecolsp][nnp];
          v_theta[ecolsp][k] = v_theta[ecolsp][nnp];
          v_z[ecolsp][k] = v_z[ecolsp][nnp];
          newvel(rengy,vel,&v_r[ecolsp][k],&v_theta[ecolsp][k],&v_z[ecolsp][k],0);
          r[ecolsp][k] = r[ecolsp][nnp];

          /* determine the charge density due to new electron */
        
          frac= (r_grid[s+1] - r[ecolsp][nnp])/dr_n[s];
          sp_n_mcc[ecolsp][s]  += frac;
          sp_n_mcc[ecolsp][s+1]+= (1.0 - frac); 

          /*-------------------------------------------------------*/
          /* create ion and assign velocities if ion species exist */

          if (ionspecies)
          {
            k = np[ionsp];
            maxwellv(&v_r[ionsp][k],&v_theta[ionsp][k],&v_z[ionsp][k],vg_therm);
            r[ionsp][k] = r[ecolsp][nnp];
  
            /* determine the charge density due to newly created ion */
            sp_n_mcc[ionsp][s]  += frac;
            sp_n_mcc[ionsp][s+1]+= (1.0 -  frac); 
          }

          /*---------------------------------------*/
          /* finally scatter the incident electron */

          vel = sqrt(engy/Escale[ecolsp]);
          newvel(engy, vel, &v_r[ecolsp][nnp], &v_theta[ecolsp][nnp], 
                 &v_z[ecolsp][nnp], 0);

          if(++np[ionsp] >maxnp[ionsp] || ++np[ecolsp] >maxnp[ecolsp])
          {	  
          #ifdef MPI_1D
	          if (my_rank == ROOT)
               printf("mcc(Ionization): too many ions created ");
	       #else
             printf("mcc(Ionization): too many ions created ");
	       #endif
             exit(1);
          }

          if(theRunWithXFlag)
          {
            mccrate[5][s]  += frac/vol[s]/icorrect;
            mccrate[5][s+1]+= (1.0 - frac)/vol[s+1]/icorrect; 
          }

          /*---------------------------------------*/
          /*      Reduce excited state density     */
          k=s/nc_ratio;
          frac=sp_ex[k]/(sp_ex[k]+sp_exm[k]);
          sp_ex[k]-=frac/vol[s];
          sp_exm[k]-=(1-frac)/vol[s];

        }   /* end of test for type of collision */
      }    /* end of for(j=0; j< N; j++)         */
      icorrect=0;
    }     /* end of if (dt_correct<=icorrect)    */
  }  
  /*------------------------------------------------------------* 
             Electron Collision : Ar  + e  
   *------------------------------------------------------------*/
  if(ecollisional && isp==ecolsp)
  {
    if(theRunWithXFlag) {
      for (j=0; j<ng; j++) {
        mccrate[0][j] = mccrate[1][j] = mccrate[2][j] = mccrate[6][j] =
        Pel[j] = Pie[j] = 0.0;
      }
     
      nu_ex = nu_el = nu_iz = 0;
   #ifdef POW
      els_loss = 0.0;
   #endif
    }

 #ifdef POW
    else {
      els_loss = 0.0;
      nu_ex = nu_iz = 0;
    }
 #endif

    ecol_extra += np[ecolsp]*col_prob_e;
    N = ecol_extra;
    ecol_extra -= N;
    
    nnp = np[ecolsp];

    for(j=0; j< N; j++)
    {
      nnp--;
      index= nnp*frand();

/* Move colliding particle to end of array, so it isn't chosen twice */
      swap(r[ecolsp],index,nnp);
      swap(v_r[ecolsp],index,nnp);
      swap(v_theta[ecolsp],index,nnp);
      swap(v_z[ecolsp],index,nnp); 

   /*------------------------------------------------------------*/
      /* determine if a collision occurs */
      v_sqrd = v_r[ecolsp][nnp]*v_r[ecolsp][nnp] +v_theta[ecolsp][nnp]*v_theta[ecolsp][nnp]
	     +v_z[ecolsp][nnp]*v_z[ecolsp][nnp];
      engy= Escale[ecolsp]*v_sqrd;
      vel = sqrt(v_sqrd);
      sigma_total = max_sigmav_e/(vel*vscale);
      rand_num= frand();
      
      /*------------------------------------------------------------*/
      /* determine the type of collision and calculate new velocities, etc */
      /*------------------------------------------------------------*/
      /* if the collision is elastic */
      
      if (rand_num <= (sum_prob=asigma1(engy)/sigma_total) )
      {
	/* first normalize the velocity components by the speed, vel */
	v_r[ecolsp][nnp] /= vel;
	v_theta[ecolsp][nnp] /= vel;
	v_z[ecolsp][nnp] /= vel;
	
	/* scatter the electron */
	newvel(engy, vel, &v_r[ecolsp][nnp], &v_theta[ecolsp][nnp], &v_z[ecolsp][nnp], 1);

        /* collect diagnostics*/
        if(theRunWithXFlag)
        {
          s = j_pointer(ecolsp,nnp);
	  frac= (r_grid[s+1] - r[ecolsp][nnp])/dr_n[s];
	  temp = frac/vol[s];
	  mccrate[0][s]  += temp;
	  Pel[s] += temp*en_loss;
	  temp = (1.0 - frac)/vol[s+1]; 
	  mccrate[0][s+1]+= temp;
	  Pel[s+1] += temp*en_loss;
          nu_el ++; 
       #ifdef POW
          els_loss += en_loss;
       #endif
        }
     #ifdef POW
        else els_loss += en_loss;
     #endif
      }
      /**********************************/
      /* Excitation collision           */ 
      /* First check for metastable collisions using asigma7 (which is a subset of the TOTAL excitation x-s asigma2). 
       Excitation to the radiative state (effective x-sect = asigma2 - asigm7) will be tested in the next section. 
       nb sum_prob is NOT increased here, but in next section, where the total excitation probability is used. This 
       format requires metastable excitation to be tested before radiative excitation */

      else if (engy>e_me && rand_num<=(sum_prob + asigma7(engy)/sigma_total))
      {
        /* first normalize the velocity components by the speed, vel */
        v_r[ecolsp][nnp] /= vel;
        v_theta[ecolsp][nnp] /= vel;
        v_z[ecolsp][nnp] /= vel;

        /* subtract the excitation energy */
        engy -= e_me;
        vel   = sqrt(fabs(engy)/Escale[ecolsp]);

        /* scatter the electron */
        newvel(engy, vel, &v_r[ecolsp][nnp], &v_theta[ecolsp][nnp], &v_z[ecolsp][nnp], 0);

        /* Determine the excitation profile */
        if(theRunWithXFlag)
        {
          s = j_pointer(ecolsp,nnp);
          frac= (r_grid[s+1] - r[ecolsp][nnp])/dr_n[s];
	  temp = frac/vol[s];
          mccrate[6][s]  += temp;
	  Pie[s]  += e_me*temp;
	  temp = (1.0 - frac)/vol[s+1]; 
          mccrate[6][s+1]+= temp;
	  Pie[s+1]+= e_me*temp; 
          nu_ex ++; 
        }
     #ifdef POW
        else nu_ex++;
     #endif 
      }

      /*------------------------------------------------------------*/
      /* Now check for radiative excitation and increase sum_prob by total
         excitation probability */
      
      else if (engy > eexc && rand_num <= (sum_prob += asigma2(engy)/sigma_total))
      {
	/* first normalize the velocity components by the speed, vel */
	v_r[ecolsp][nnp] /= vel;
	v_theta[ecolsp][nnp] /= vel;
	v_z[ecolsp][nnp] /= vel;
	
	/* subtract the excitation energy */
	engy -= eexc;
	vel   = sqrt(fabs(engy)/Escale[ecolsp]);
	
	/* scatter the electron */
	newvel(engy, vel, &v_r[ecolsp][nnp], &v_theta[ecolsp][nnp], &v_z[ecolsp][nnp], 0);
	
	/* Determine the excitation profile */