diagnostics.c

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

/**********************************************************************
** This file contains time-history and frequency diagnostics. Includes:
** TIME HISTORY
**   mccrate - collision rates
**   
***********************************************************************/
#include "pdc1.h"

/************************************************************************/
/* time history accumulator; calculates and stores all history values, and
   performs combing on history values when low on memory */

FILE *Output;

void freqanalysis(void);
void realft(float [], int, int);
void velocity(void);
float mu_calc(void);

int sheath0(void)
{
  int ii;
  for (ii = 0; ii < nc/2; ii++) {
    if (sp_n[1][ii] && sp_n[1][ii+1])
     if (sp_n[0][ii]/sp_n[1][ii] + sp_n[0][ii+1]/sp_n[1][ii+1] >= 1.2) break;
  }
  return ii;
}

int sheath1(void)
{
  int ii;
  for (ii = nc; ii > nc/2; ii--) {
    if (sp_n[1][ii] && sp_n[1][ii-1])
      if (sp_n[0][ii]/sp_n[1][ii] + sp_n[0][ii-1]/sp_n[1][ii-1] >= 1.2) break;
  }
  return nc-ii;
}

void history()
{
  int i, ii, j, k,isp;
  static float sntos, fscale, ef_norm;
  float power, temp, temq, Iz_2;
  float vz_sum;

  static int count=1, navecount=0;
  static int d_skip=1;
  /* static int icount=0; */
  static float pow_avg=0.0, jnorm[NSMAX];
  static float mu_3;

  /* Calculate normalisation constants */
  if (rescale_on) 
  {
    sntos=ECHARGE*nc2p/(2.0*PI*dr*height); /* Normalisation for sigma diag */
    fscale = 0.5*PI*height*dr*dr/ptopn/ECHARGE;
    for (isp=0; isp<nsp; isp++)
      jnorm[isp] = charge[isp]*nc2p/(PI*dr*height*dt);
    ef_norm = 1.0/(ptopn*dr);
    rescale_on = 0;
  }

/* Renormalise current */
  I_rf*=intoi[0];

/* Calculate power  */
  power = I_rf*phi[0]/ptopn;                
  if (it[0]) pow_avg = (pow_avg*(it[0]-1) + power)/((float)it[0]);

/* Normalise charge density */
  for(i=0; i<ng; i++) 
  {
	 rho[i] *= dntod*ECHARGE;
	 rho_unfilt[i] *= dntod*ECHARGE;
  }

/* Normalise density for plotting and calculate ke and j profiles*/
  for(isp=0; isp<nsp; isp++)
  {
    for(i=0; i<ng; i++) density[isp][i] = sp_n[isp][i]*dntod;
    if (k_count[isp] == 1)
    {
      for(i=0; i<ng; i++) 
      {
         ke_r[isp][i] *= Escale[isp];
/*
         ker_r[isp][i] *= Escale[isp];
         ket_r[isp][i] *= Escale[isp];
         kez_r[isp][i] *= Escale[isp];
*/
         j_r[isp][i] *= jnorm[isp];
         j_z[isp][i] *= jnorm[isp];
/*
         jr_unfilt[isp][i] *= jnorm[isp];
         jz_unfilt[isp][i] *= jnorm[isp];
*/
      }

/**** Store angular momentum */
      for (i=0; i<np[isp]; i++) 
         angular_mtm[isp][i] = v_theta[isp][i]*r[isp][i];
    } /* end if (k_count ...) loop */
  }

/* Calc total axial current using both sum of v_z (1) and integration of J_z (2) */
  Iz_2 = 0.0;
  for (isp=0; isp<nsp; isp++)
      for (ii=0; ii<ng; ii++) Iz_2 += j_z[isp][ii]*r_array[ii]*dr_n[ii];
  Iz_2 *= TWOPI*dr;
  
  /*if (src !='P' && src !='p' && src !='J' && src!='j')  
  {
      vz_sum = 0;
      for (ii=0; ii<np[0]; ii++) vz_sum += v_z[0][ii];
      I_z = charge[0]*vscale*vz_sum*nc2p/height;
  }*/

/********* Calculate time-averaged profiles ******************/
  if (RT_flag && it_rate) {
    for(i=0; i<2; i++)
      for (j=0; j<nc_RT; j++)
        ex_rate_show[i][j] = (ex_rate_show[i][j]*(it_rate-1)
                           + ex_mccrate[i][j])/((float)it_rate);
    for(i=0;i<nc_RT;i++) {
      prod[i]+=ex_rate_show[0][i];
      prodm[i]+=ex_rate_show[1][i];
    }
  }
  if (it[0])
  {
/* Calculate collision profiles  */
     if (ecollisional || icollisional)
     {	
       for (i=0; i<ncolls; i++)
         for (j=0; j<ng; j++) 
           rate_show[i][j] = (rate_show[i][j]*(it[0]-1) +
                             mccrate[i][j])/((float)it[0]);
          /* c.f CVS version,
           rate_show[i][j] = (rate_show[i][j]*(it[0]-1) +
                        mccrate[i][j]*dntod/dt)/((float)it[0]);
          */
       for (j=0; j<ng; j++) {
         Pel_show[j] = (Pel_show[j]*(it[0]-1) + 
                       ECHARGE*Pel[j]*dntod/dt)/((float)it[0]);
         Pie_show[j] = (Pie_show[j]*(it[0]-1) +
                       ECHARGE*Pie[j]*dntod/dt)/((float)it[0]);
         Pcx_show[j] = (Pcx_show[j]*(it[0]-1) +
                       ECHARGE*Pcx[j]*dntod/dt)/((float)it[0]);
       }
  

/*       if (RT_flag) for (i=0;i<nc_RT;i++) 
       {
         k=i*nc_ratio;
         temp=temq=0;
         for (j=0;j<nc_ratio;j++) 
         {
           temp+=rate_show[1][k]+rate_show[1][k+1];   
           temq+=rate_show[6][k]+rate_show[6][k+1]; 
           k++;
         }
         prod[i]+=0.5*temp/nc_ratio;
         prodm[i]+=0.5*temq/nc_ratio;
       }   */
     }

/* Calculate time-averaged profiles - J.E, K.E., J */
    for (isp=0;isp<nsp;isp++)
    {
      if (k_count[isp] == 1)
      {
/* Radial & Axial Current density profile */
         for (j=0; j<ng; j++) 
         {
            jr_prof[isp][j] = (jr_prof[isp][j]*(it[isp]-1) + j_r[isp][j])/it[isp]; 
            jz_prof[isp][j] = (jz_prof[isp][j]*(it[isp]-1) + j_z[isp][j])/it[isp]; 
         }

/* Calculate electron & ion mobility profiles from Ez and Jz */
       for (i=0; i<ng; i++)
	 {
	   if(density[isp][i])
	     {
	       if(E_z > 1E-5)
		 {
		   mu_prof[isp][i]=(mu_prof[isp][i]*(it[isp]-1) + jz_prof[isp][i]/(ECHARGE*density[isp][i]*E_z*ef_norm))/it[isp];
		 }
	       else
		 {
		   mu_prof[isp][i] = 0.0;
		 }
	     }
	 }

/* Time-averaged J.E */
         for (j=0; j<ng; j++) 
         {
	     jdote[isp][j] = (jdote[isp][j]*(it[isp]-1) + j_r[isp][j]*efield[j]*ef_norm)/((float)it[isp]);
             jzdotez[isp][j] = (jzdotez[isp][j]*(it[isp]-1) + j_z[isp][j]*E_z*ef_norm)/((float)it[isp]);
	 }
      } /* end of if k_count = 1 */

/* Average kinetic energy profile */
         for (j=0; j<ng; j++)
         {
              ke_avg[isp][j] = (ke_avg[isp][j]*(it[0]-1) + ke_r[isp][j])/it[0];
/*
              ker_avg[isp][j] = (ker_avg[isp][j]*(it[0]-1) + ker_r[isp][j])/it[0];
              ket_avg[isp][j] = (ket_avg[isp][j]*(it[0]-1) + ket_r[isp][j])/it[0];
              kez_avg[isp][j] = (kez_avg[isp][j]*(it[0]-1) + kez_r[isp][j])/it[0];
*/
         }


/* Time-averaged mid-plasma energy distribution function */
       for (i=0; i<nbin_mid[isp]; i++)
          fe_mid_show[isp][i] = (fe_mid_show[isp][i]*(it[0]-1) +
                   fe_mid[isp][i]/sqrt(emid_array[isp][i]+DBL_MIN))/((float)it[0]);
    }
  }

/* Determine velocity distribution */
  if (vel_dist_accum) velocity();
  
/* Determine electron mobility using 3rd method */
  /*if (k_count[0] == 1) mu_3 = mu_calc();*/

/* Store evolution of current, mid and electrode potential and power for fft */
  if(nfft)
  {
      if (thist_hi >= nfft)
      {
	  freqanalysis();
	  thist_hi=0;
      }
      cur_hist[thist_hi]= I_rf;
      phi_hist[0][thist_hi] = phi[0];
      phi_hist[1][thist_hi] = phi[nc/2];
      pow_hist[thist_hi] = power;
      Local_t_array[thist_hi] = atime;
      thist_hi++;

  /*- Time-averaged species densities -----*/

      if(navecount == nfft)     
      {
        for(isp=0; isp<nsp; isp++)
          for(i=0; i<ng; i++)
          {
            den_ave[isp][i]= den_sum[isp][i]/nfft; 
            den_sum[isp][i]= density[isp][i];
          }

        for(i=0; i<ng; i++)
        {
  /*- Average charge density profile -----*/
           rho_ave[i] = rho_sum[i]/nfft;
           rho_sum[i] = rho[i];

  /*- Time-averaged potential and electric field -----*/
          e_ave[i]= e_sum[i]/nfft;      
          e_sum[i]  = efield[i];
          phi_ave[i]= phi_sum[i]/nfft;  
          phi_sum[i]= phi[i];
        }
        navecount =0;
     }
 /* Sum up quantities used for doing averages */
    else
    {

/* Sum density for plotting average */
      for(isp=0; isp<nsp; isp++)
        for(i=0; i<ng; i++) den_sum[isp][i]+= density[isp][i];

/* Sum charge density */
      for(i=0; i<ng; i++) rho_sum[i]+= rho[i];

 /*- Instantaneous potential and E field profiles ----*/
      for (i=0; i<ng; i++)
      {
        e_sum[i] += efield[i];
        phi_sum[i] += phi[i];
      }
      navecount++;
    }    /* end of "else  */
  }      /* end of "if (nfft)" */

  /*------------Check # diagnostics stored -------------------------*/

  if (--count) return;	  /* only accum every d_skip steps */
  if (hist_hi >= HISTMAX) /* comb time histories */
  {
     for (isp=0; isp<nsp; isp++)
     {
	for (i=1, k=4; i<HISTMAX/4; i++, k+=4)
	{
           np_hist[isp][i] = np_hist[isp][k];
	   np_hist2[isp][i] = np_hist[isp][k];
           kes_hist[isp][i] = kes_hist[isp][k];
	}
     }
     for (i=1, k=4; i<HISTMAX/4; i++, k+=4)
     {
	  com_phi_hist[0][i] = com_phi_hist[0][k];
	  com_phi_hist[1][i] = com_phi_hist[1][k];
	  com_power_hist[i] = com_power_hist[k];
	  avg_pow_hist[i] = avg_pow_hist[k];
	  com_cur_hist[i] = com_cur_hist[k];
	  sigma_hist[i] = sigma_hist[k];
	  ese_hist[i] = ese_hist[k];
          sheath[0][i] = sheath[0][k];
          sheath[1][i] = sheath[1][k];
	  ke_tot[i] = ke_tot[k];
	  te_hist[i] = te_hist[k];
          Ez_hist[i] = Ez_hist[k];
          Iz_hist[i] = Iz_hist[k];
          Iz_hist2[i] = Iz_hist2[k];
          mu_hist[i] = mu_hist[k];
          mu_hist2[i] = mu_hist2[k];
/*          mu_hist3[i] = mu_hist3[k];*/
/*          nu_iz_hist[i] = nu_iz_hist[k]; */
/*          nu_ex_hist[i] = nu_ex_hist[k]; */
/*          nu_el_hist[i] = nu_el_hist[k]; */
	  t_array[i] = t_array[k];
          C_hist[0][i] = C_hist[0][k];			/* LEEHJ */
          C_hist[1][i] = C_hist[1][k];
     }
     hist_hi = i;