rescale.c

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

/***********************************************************************
** This routine rescales the particle numbers when they get too large
** or too small. For reducing the number of simulation particles every
** n th particle is removed, and then nc2p is adjusted so that the number
** of real particles represented is constant. When particles are added
** every mth particle is "doubled" and added to the end of the array.
**
** After the particle numbers are adjusted and nc2p is rescaled parameters
** which are affected by nc2p must be scaled also. ie charge and current
** which are normalised by nc2p, the alphas used for calculating the potential,
** and some of the normalisation constants.
**  
**      in:	npart
***********************************************************************/
#include "pdc1.h"

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

void rescale(int k)
{ 
  int ii, jj, isp;
  int nnp, np_old;
  float nc2p_frac;

/* If the number of simulation particles is too large */
  np_old = np[k];
  if (np_old > 0.9*maxnp[0]) 
  {
  #ifdef MPI_1D
    if(my_rank == ROOT){
      printf("Particle number too large, rescaling nc2p \n");
      printf("     old npart = %d %d, old nc2p = %e \n ",np[0], np[1], nc2p);
      printf("     t = %e it = %d ",atime, it[0]);
    }
  #else
    printf("Particle number too large, rescaling nc2p \n");
    printf("     old npart = %d %d, old nc2p = %e \n ",np[0], np[1], nc2p);
    printf("     t = %e it = %d ",atime, it[0]);
  #endif
/* Sort particles by position */
    for (jj=0; jj<nsp; jj++) 
         sort(np[jj], r[jj]-1, v_r[jj]-1, v_theta[jj]-1, v_z[jj]-1);

    for (isp=0; isp<nsp; isp++)
    {
        ii = 1;
        nnp = 2*((np[isp]-1)/2);
        do
        {
           r[isp][ii] = r[isp][nnp];
           v_r[isp][ii] = v_r[isp][nnp];
           v_theta[isp][ii] = v_theta[isp][nnp];
           v_z[isp][ii] = v_z[isp][nnp];
           nnp -= 2;
           ii += 2;
        }  while(ii<nnp);
        np[isp]=nnp+1;
    }

/* If the number of simulation particles is too small */
  } 
  else if (np_old < min_np)
  {
  #ifdef MPI_1D
    if(my_rank == ROOT){
      printf("Particle number too small, rescaling nc2p \n");
      printf("     old npart = %d %d, old nc2p = %e \n ",np[0], np[1], nc2p);
    }
  #else
    printf("Particle number too small, rescaling nc2p \n");
    printf("     old npart = %d %d, old nc2p = %e \n ",np[0], np[1], nc2p);
  #endif

    /* loop over each species */
    for (isp=0; isp<nsp; isp++)
    {
       for (ii=0; ii<np[isp]; ii++)
       {
          r[isp][ii+np[isp]] = r[isp][ii];
          v_r[isp][ii+np[isp]]= v_r[isp][ii];
          v_theta[isp][ii+np[isp]]= v_theta[isp][ii];
          v_z[isp][ii+np[isp]]= v_z[isp][ii];
    }

        np[isp]*=2;
    }
  }
  nc2p_frac = (float)np_old/np[k];

  /* Rescale parameters which contain nc2p explicitly or implicitly */
  nc2p = nc2p*nc2p_frac;

#ifdef MPI_1D
  if(my_rank == ROOT)
    printf("     new npart = %d %d new value nc2p = %e \n ",np[0], np[1], nc2p); 
#else
    printf("     new npart = %d %d new value nc2p = %e \n ",np[0], np[1], nc2p); 
#endif

  rhon *= nc2p_frac;  
  for (isp=0; isp<nsp; isp++) intoi[isp] *= nc2p_frac;
  alpha0 *= nc2p_frac;
  alpha1 *= nc2p_frac;
  alpha2 *= nc2p_frac;
  alpha3 *= nc2p_frac;
  alpha4 *= nc2p_frac;
  sigma /= nc2p_frac;
  oldsigma /= nc2p_frac;  
  q_0 /= nc2p_frac;    
  q_1 /= nc2p_frac;    
  q_2 /= nc2p_frac;    
  q_3 /= nc2p_frac;    
  I_rf /= nc2p_frac;
  if(src== 'I' || src== 'i' || src == 'K')
  {
      acbias /=nc2p_frac;
      dcbias /=nc2p_frac;
      acbias2 /=nc2p_frac;
      dcbias2 /=nc2p_frac;
  } 
  if (src=='J' || src=='j')
  { 
      acbias /=nc2p_frac;
      dcbias /=nc2p_frac;
      E_z /= nc2p_frac;
  }


  rescale_on=1;    /* Turn rescale indicator on */

} /* end rescale */
/****************************************************************/
/*  Heap sort  */
void sort(int n, float *ra, float *rb, float *rc, float *rd)
{
  int l, j, ir, i;
  float rra, rrb, rrc, rrd;

  l = (n >>1) +1;
  ir= n;

  for(;;)
  {
    if(l >1)
    {
      rra = ra[--l];
      rrb = rb[l];
      rrc = rc[l];
      rrd = rd[l];
    }
    else
    {
      rra   = ra[ir];
      rrb   = rb[ir];
      rrc   = rc[ir];
      rrd   = rd[ir];

      ra[ir]= ra[1];
      rb[ir]= rb[1];
      rc[ir]= rc[1];
      rd[ir]= rd[1];


      if(--ir == 1)
      {
        ra[1] = rra;
        rb[1] = rrb;
        rc[1] = rrc;
        rd[1] = rrd;
        return;
      }
    }
    i=l;
    j= l << 1;
    while(j <= ir)
    {
      if(j < ir && ra[j] < ra[j+1]) ++j;
      if(rra <ra[j])
      {
        ra[i] = ra[j];
        rb[i] = rb[j];
        rc[i] = rc[j];
        rd[i] = rd[j];
        j += (i=j);
      }
      else j=ir+1;
    }
    ra[i]= rra;
    rb[i]= rrb;
    rc[i]= rrc;
    rd[i]= rrd;
  }
}
/****************************************************************/