field.c

XPDC1是针对静电模型中的等离子体模拟程序！通过INP文件输入参数有较好的通用性

/*******************************************************************/
/*	This set of routines does the charge weighting to the      */
/*	grid, calculates the potentials and fields at the grid,    */
/*	interpolates the fields back to the charged species and    */
/*	then moves them. 					   */
/*                                                                 */
/*	* the weighting can be linear or volume                    */
/*      * potentials are solved using a backwards tridiagonal      */
/*	matrix method.    					   */
/*                                                                 */
/*      It is based on the routine pflmv.c, but has been           */
/*      re-written                                                 */
/*      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) The charge density is smoothed using a one-to-one       */
/*         filter to noise in the fields                           */
/******************************************************************/
#include "pdc1.h"
#include <xgrafix.h>

float source(void);
void circuit1(void), circuit2(void), circuit3(void), circuit4(void),
circuit5(void), circuit6(void), circuit7(void); 
void  (*circuitptr)(void);
void DigitalFilter(float *, float *, int, int);

float *xp, *a_phi, *b_phi, *c_phi;
float ke_sum, ke_old; 
float Iz_norm, mu_norm;
float vz_sum, sv_sum;
#ifdef MPI_1D
float sum[2], sum_tot[2];
#endif
/***************************************************************/
void gather(int isp)
/* Accumulate charge density on grid */
{
  int ii, j, s;
  float frac;
  float vr_sq, vt_sq, vz_sq, energy, vel, sv_avg;
  
/* Set current, ke profile and fe diagnostics to 0 */
    if (theRunWithXFlag)
    {
      for (ii=0; ii<=nc; ii++) {
         j_r[isp][ii]=j_z[isp][ii]=ke_r[isp][ii]=0.0;
         ker_r[isp][ii]=ket_r[isp][ii]=kez_r[isp][ii]=0.0; 
      }
      ke_sp[isp] = 0.0;
      for (ii=0; ii<=nbin_mid[isp]; ii++) fe_mid[isp][ii]=0.0;
      sv_avg = 0;
    }

/* Calculate total density at t-1 due to position + ioniz colls and 0
sp_n_mcc and sp_n_k */
  for (ii=0; ii<=nc; ii++) 
  {
     sp_n_0[isp][ii] = sp_n_k[isp][ii] + sp_n_mcc[isp][ii];
     sp_n_mcc[isp][ii]= sp_n_k[isp][ii]= 0.0;
  }

/* Loop over number of particles */
    for (ii=0; ii<np[isp]; ii++) 
    {
/* check that particle is between electrodes */
      if (r[isp][ii]>=r0 && r[isp][ii]<r1) 
      {
        j = j_pointer(isp,ii);

        /* linear weighting */
	frac = (r_grid[j+1] - r[isp][ii])/dr_n[j];
	sp_n_k[isp][j] += frac;
	sp_n_k[isp][j+1] += (1.0 - frac);

/*   Calculate total ke of species, store kinetic energy and current at grid
cells for profile diagnostics and determine electron energy distribution in
plasma  */
        if (theRunWithXFlag)
        {
           j_r[isp][j] += frac*v_r[isp][ii];
           j_r[isp][j+1] += (1.0 - frac)*v_r[isp][ii];
           j_z[isp][j] += frac*v_z[isp][ii];
           j_z[isp][j+1] += (1.0 - frac)*v_z[isp][ii];
           vr_sq = v_r[isp][ii]*v_r[isp][ii];
           vt_sq = v_theta[isp][ii]*v_theta[isp][ii];
           vz_sq = v_z[isp][ii]*v_z[isp][ii];
           energy = vr_sq + vt_sq + vz_sq;
           ke_sp[isp] += energy;
           ke_r[isp][j] += frac*energy;
           ke_r[isp][j+1] += (1.0 - frac)*energy;
           ker_r[isp][j] += frac*vr_sq;
           ker_r[isp][j+1] += (1.0 - frac)*vr_sq;
           ket_r[isp][j] += frac*vt_sq;
           ket_r[isp][j+1] += (1.0 - frac)*vt_sq;
           kez_r[isp][j] += frac*vz_sq;
           kez_r[isp][j+1] += (1.0 - frac)*vz_sq;
	   
           if ((r[isp][ii] - rs_mid[isp])*(rf_mid[isp] - r[isp][ii]) >= 0)
           {
              s = (energy - emin_mid[isp])/de_mid[isp];
              if (0 <= s && s < nbin_mid[isp]) fe_mid[isp][s] += nc2p;
           }

/* Calculate average coll freq if NOT using axial current as source */
           if (isp==0 && src != 'J' && src != 'j')
           {
             vel = sqrt(energy);
             energy *= Escale[0];
             sv_avg += (*sigmaTptr)(energy)*vel;
           }
        } /* end of diagnostic loop */
      } /* end of "if (r[isp][i] > r0" loop */
    }   /* end of  loop over all particles */
    
/* Divide ke_r by particle number to get average ke per particle */
    if(theRunWithXFlag) {
      for (ii=0; ii<ng; ii++) {
        if (sp_n_k[isp][ii]) 
        {
           ke_r[isp][ii] /= sp_n_k[isp][ii];
           ker_r[isp][ii] /= sp_n_k[isp][ii];
           ket_r[isp][ii] /= sp_n_k[isp][ii];
           kez_r[isp][ii] /= sp_n_k[isp][ii];
        }
	
	j_r[isp][ii] /= vol[ii];
	j_z[isp][ii] /= vol[ii];
      }
      
     /* Calculate theoretical electon mobility */
     if (isp==0 && src !='J' && src !='j') mu_e = mu_norm*np[0]/sv_avg;

     } 
}
/***************************************************************/
/*  The init_rho rouutine is used to determine the charge density */
/*  at the initial (t=0) step                                   */
/***************************************************************/
void init_rho(void)
{
  int j, isp;

  for(isp=0; isp<nsp; isp++)
  {
    it[isp] = 0;
    k_count[isp] = 0;
    gather(isp);
    for (j=0; j<ng; j++)
    {
      sp_n[isp][j]= sp_n_k[isp][j];
      sp_n_mcc[isp][j]= 0.0;
    }
  }
}
/***************************************************************/
void fields()
{
  int isp, ii;

/* Find rho(x) = total charge density...
     ..first put in background (rigid) charge density */
  for (ii=0; ii<=nc; ii++) rho[ii]= rhoback;

/* ... now add in species charge densities */
  for (isp=0; isp<nsp; isp++)
    for (ii=0; ii<=nc; ii++) 
      #ifdef MPI_1D	
        rho[ii] += q_ratio[isp]*sp_n_tot[isp][ii];
      #else
        rho[ii] += q_ratio[isp]*sp_n[isp][ii];
      #endif

/* Digitally smooth the net charge density */
  if (nsmoothing > 0) 
  {
	 for (ii = 0; ii<ng; ii++) rho_unfilt[ii] = rho[ii];
	 DigitalFilter(rho_unfilt, rho, ng, nsmoothing);
  }

  /* determine convective (particle) current to powered electrode */
  i_conv = backgrnd_i;                                /* first background current */
  #ifdef MPI_1D
    for (isp=0; isp< nsp; isp++) i_conv += iwall_tot[isp];  /* now species current */
   #else
    for (isp=0; isp< nsp; isp++) i_conv += iwall[isp];  /* now species current */
   #endif 

  /* Solve a tridiagnal matrix to get potential at time t*/
  (*circuitptr)();  /* pointer to appropriate circuit calc */ 
  
  /* Now calculate electric field at grid from potentials */
  for (ii=1; ii<nc; ii++) efield[ii] = (phi[ii-1] - phi[ii+1])/(dr_n[ii-1] + dr_n[ii]); 
  if (r0 > 0) efield[0] = ((r0+0.5*dr_n[0])*(phi[0]-phi[1]) - 0.5*rhon*rho[0]*vol[0])/r0;
  else efield[0] = 0.0;
  efield[nc]=((r1-0.5*dr_n[nc-1])*(phi[nc-1]-phi[nc])/dr_n[nc-1] + 0.5*rhon*rho[nc]*vol[nc])/r1;

}	/* end of FIELDS */

/**************************************************************/
/* Determines time-dependent source voltage. Can have         */
/*    1) dc voltage with a linear ramp                        */
/*    2) ac voltage with an exponential ramp                  */
/*    3) a combination of both                                */
/* It doesn't do a cosine ramp for the dc voltage anymore     */ 
/**************************************************************/
float source(void)
{
  float Sdc, Sac;
  float Sdc2, Sac2;

  Sdc = dcbias;
  Sac = acbias*sin(atime*w0 + theta0);

  if (ramped)
  {
    if(dcbias && atime < risetime) Sdc = (dc_ramp*atime);
    if(ac_ramp > 0.0) Sac *= exp(-ac_ramp/atime);
    if(ac_ramp < 0.0) Sac *= exp(atime/ac_ramp);
  }
  if((src == 'W') ||(src == 'K'))
    {
      Sdc2 = dcbias2;
      Sac2 = acbias2*sin(atime*w02 + theta02);
      
      if (ramped2)
        {
          if(dcbias2 && atime < risetime2) Sdc = (dc_ramp2*atime);
          if(ac_ramp2 > 0.0) Sac *= exp(-ac_ramp2/atime);
          if(ac_ramp2 < 0.0) Sac *= exp(atime/ac_ramp2);
        }
      
      return (Sac2 + Sdc2 + Sac + Sdc);
      
    }

  return (Sac + Sdc);
}

/***************************************************************/
/* Smoothing the array using a digital filter analogous to 1-2-1
** method used in planar system with the proper boundary conditions 
** to conserve charge. Note that nSmoothing is only used on first
** pass thru to set up coefficients for multipass smoothing */

void DigitalFilter(float *source, float *dest, int N, int nSmoothing)
{
  static int init = 0;
  static float *alpha, **A;
  int j, c, d;
  float Ajd, Ajdm1;
  if (nSmoothing < 1) return;
  N--; /* last element index: 0...N */

  /** Initialize arrays on first use of filter **/
  if (init==0)
  {
	 int i;
	 init = 1;
	 alpha = (float *) malloc((N+1)*sizeof(float));
	 alpha[0] = 0.5;
	 alpha[1] = alpha[0]*vol[0]/vol[1];
	 for (j=1; j<N; j++)
		alpha[j+1] = (2*alpha[j]*vol[j] - alpha[j-1]*vol[j-1])/vol[j+1];
	 A = (float **) malloc((N+1)*sizeof (float *));
	 for (j=0; j<=N; j++) {
		A[j] = (float *) calloc((2*nSmoothing + 1),sizeof(float));
		/* intitial values for one pass */
		for (c=max(0,j-1); c<=min(N,j+1); c++){
		  d = c - j + nSmoothing;
		  if (j==c) 
			 if (j==0 || j==N) A[j][d] = 1 - alpha[j];
			 else A[j][d] = 1 - 2*alpha[j];
		  else A[j][d] = alpha[j];
		}
	 }
	 for (i=1; i<nSmoothing; i++) {
		for (j=0; j<=N; j++) {
		  Ajd = A[j][0];
		  for (c=max(0,j-nSmoothing); c<=min(N,j+nSmoothing); c++){
			 d = c - j + nSmoothing;
			 Ajdm1 = Ajd; /* shift */
			 Ajd = A[j][d]; /* save new element */
			 if (c == 0) A[j][d] = Ajd*(1-alpha[c]) + A[j][d+1]*alpha[c+1];
			 else if (c==N) A[j][d] = Ajdm1*alpha[c-1] + Ajd*(1-alpha[c]);
			 else A[j][d] = Ajdm1*alpha[c-1] + Ajd*(1-2*alpha[c]) + A[j][d+1]*alpha[c+1];
		  }
		}
	 }
  }
  for (j=0; j<=N; j++) {
	 dest[j] = 0;
	 for (c=max(0,j-nSmoothing); c<=min(N,j+nSmoothing); c++){
		d = c - j + nSmoothing;
		dest[j] += source[c]*A[j][d];
	 }
  }
}

/***************************************************************/
/*    Initialise parameters for calculating the field solve    */
/*    using tridiagonal matrix method                          */
/***************************************************************/
void field_init(void)
{
  int i, k;
  float indn, capn, resn;
  float jz_norm, pz_norm, r_avg;

/****** Create space for alphas, as, cs, etc */
  a_phi = (float  *)malloc(nc*sizeof(float));
  b_phi = (float  *)malloc(nc*sizeof(float));
  c_phi = (float  *)malloc(nc*sizeof(float));
  xp = (float  *)malloc(nc*sizeof(float));

  if(!a_phi || !b_phi || !c_phi || !xp )
  {
  #ifdef MPI_1D
    if(my_rank == ROOT)
      printf("start() : malloc() a, c, xp returned NULL");
  #else
    printf("start() : malloc() a, c, xp returned NULL");
  #endif
    exit(-1);
  }

/* normalise circuit parameters */
     capn = extc/(ptopn*ECHARGE*nc2p);
     indn = extl*ptopn*ECHARGE*nc2p/dt/dt;
     resn = extr*ptopn*ECHARGE*nc2p/dt;

/**Calc constants for pot calc */
  for (i=1; i<nc; i++)
  {
     r_avg = (dr_n[i] + dr_n[i-1])*0.5;
     a_phi[i]=1.0/dr_n[i-1]/r_avg - dr_n[i]/(r_grid[i]*2.0*r_avg*dr_n[i-1]);
     b_phi[i]=2.0/dr_n[i]/dr_n[i-1] - (dr_n[i] - dr_n[i-1])/r_grid[i]/dr_n[i]/dr_n[i-1];
     c_phi[i]=1.0/r_avg/dr_n[i] + dr_n[i-1]/(r_grid[i]*2.0*r_avg*dr_n[i]);