field.c

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

  }
  c_phi[0] = 1.0;

/* Calc mu_norm if diagnostics turned on */
  if(theRunWithXFlag)
    mu_norm =  (*mueFacPtr)()*ECHARGE/(mass[0]*gas_den*vscale);

 /*---------------------------------------------*/
 /* Normalise external circuit parameters       */
  if(src== 'I' || src== 'i' || src=='K')
  {
      acbias *=dt/ECHARGE/nc2p;
      dcbias *=dt/ECHARGE/nc2p;
      acbias2 *=dt/ECHARGE/nc2p;
      dcbias2 *=dt/ECHARGE/nc2p;
  } else if (src=='V' || src=='v' || src=='W')
  {
      acbias *=ptopn;
      dcbias *=ptopn;
      acbias2 *=ptopn;
      dcbias2 *=ptopn;
  }
  else if (src=='J' || src=='j')
  {
    jz_norm = height*dt/(nc2p*mass[0]*vscale);
    mu_norm = (*mueFacPtr)()*ECHARGE/(mass[0]*gas_den*vscale);
    dcbias *= jz_norm;
    acbias *= jz_norm;
    circuitptr= circuit5;
    Iz_norm = charge[0]*nc2p*vscale/height;
    
    ke_sum = 0;
    for (i=0; i<np[0]; i++)
        ke_sum += v_r[0][i]*v_r[0][i] + v_theta[0][i]*v_theta[0][i] + v_z[0][i]*v_z[0][i];
     ke_old = ke_sum;
  }
  else if (src=='P' || src=='p')
  {
/* Convert to normalised power */
     pz_norm = ptopn*dr/height;
     acbias *= pz_norm;
     dcbias *= pz_norm;
     circuitptr= circuit7;
     Iz_norm = charge[0]*nc2p*vscale/height;

     ke_sum = 0;
     for (i=0; i<np[0]; i++)
        ke_sum += v_r[0][i]*v_r[0][i] + v_theta[0][i]*v_theta[0][i] + v_z[0][i]*v_z[0][i];
     ke_old = ke_sum;
  }
  else if (src=='E' || src=='e')
  { 
    dcbias *= ptopn*dr;
    acbias *= ptopn*dr;
    circuitptr= circuit6;
  }

  if (r0 == 0.0)  /* no inner conductor */
  {
    xp[0] = 1.0;
    k = 1;
  }
  else if(extc < 1e-30)  /* Open Circuit, when C tends to 0.0 */ 
  {
    if(fabs(dcbias)+fabs(acbias) >0.0)
    {
      #ifdef MPI_1D
	if (my_rank == ROOT)
	  printf("START:The active external source is ignored since C<1e-30\n");
      #else
	printf("START:The active external source is ignored since C<1e-30\n");
      #endif
       
    }
	circuitptr = circuit1;
	k=1;
	xp[0]=1.0;
  } 
  else if(src== 'I' || src== 'i' || src =='K')       
/* When current source is applied */
  {
    circuitptr= circuit2;
    k=1;
    xp[0]= 1.;
  }
  else if(indn <1e-30 && resn <1e-30 &&        /* Short Circuit */
         (extc*log(r1/r0)/(TWOPI*height*epsilon)) >= 1e5) 
  {
    circuitptr = circuit3;
    k=2;
    xp[1]=0.5;
  }
  else  /* The general case with external voltage source */
  { 
     circuitptr = circuit4;

     alpha0 = 1.5*resn + 2.25*indn + 1.0/capn;
     alpha1 = 2.0*resn + 6.0*indn;
     alpha2 = -0.5*resn - 5.5*indn;
     alpha3 = 2.0*indn;
     alpha4 = -0.25*indn;

     k=1; 
     xp[0] = 1.0/(1.0 + 0.5*rhon/(r0 + 0.5)/alpha0);
  }
  for (i=k;i<nc;i++) xp[i] = 1.0/(b_phi[i] - a_phi[i]*c_phi[i-1]*xp[i-1]);

}       /* end field_init */

/***************************************************************/
/*  1: Open circuit case */
void circuit1()

{
  int ii;
  float yp[nc];

  oldsigma = sigma;
  sigma += i_conv/r0;

  yp[0] = xp[0]*0.5*rhon*(sigma*r0 + rho[0]*(r0 + 0.25))/(r0+0.5);
  for (ii=1; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>=0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];

}

/***************************************************************/
/* 2: Current source */

void circuit2()
{
  int ii;
  float temp;
  float yp[nc];
  
  /* Calculate charge density on electrode - second order method  */
  I_rf = source();
  temp= sigma;
  sigma = (4.*sigma -oldsigma +2.0*(i_conv+I_rf)/r0)/3.0;
  oldsigma= temp;

  /* Calculate charge density on electrode */
  yp[0] = 0.5*rhon*xp[0]*(sigma*r0 + rho[0]*(r0 + 0.25))/(r0 + 0.5);
  for (ii=1; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>=0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];

}

/***************************************************************/
/* 3: Short circuit case */

void circuit3()
{
  int ii;
  float yp[nc];

  phi[0] = source();
  yp[1] = xp[1]*(a_phi[1]*phi[0] + rhon*rho[1]);
  for (ii=2; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];

  oldsigma = sigma;
  sigma = ((r0+0.5)*(phi[0]-phi[1])*2.0/rhon - (r0+0.25)*rho[0])/r0;
  I_rf = (sigma-oldsigma)*r0 - i_conv;
}

/***************************************************************/
/* 4: General voltage source case w RLC series circuit */

void circuit4()
{
  int ii;
  float k_dash;
  float yp[nc];
  
  k_dash = alpha1*q_0 + alpha2*q_1 +alpha3*q_2 +alpha4*q_3;

  yp[0] = xp[0]*0.5*rhon*(r0*sigma + i_conv + (k_dash + source())/alpha0 - q_0
           + rho[0]*(r0 + 0.25))/(r0 + 0.5);
  for (ii=1; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>=0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];

  q_3= q_2;
  q_2= q_1;
  q_1= q_0;
  q_0 = (k_dash + source() -phi[0])/alpha0;
  
  I_rf = q_0 - q_1;
  oldsigma = sigma;
  sigma += (i_conv + I_rf)/r0;
}

/***************************************************************/
/* 5: Case with no inner conductor and axial current supplied 
** This version calculates the electron mobility from the electron momentum
** transfer collision frequency.
***************************************************************/

void circuit5()
{
  int ii;
  static int skip = 1;
  float yp[nc];
  float energy, vel, Ne;

  /* calculate average collision frequency every skip timesteps */
  if (it[0]%skip == 0)
  {  
  #ifdef MPI_1D
     if (my_rank == ROOT) {
       if ( fabs(ke_old - ke_sum)/ke_old < 0.0001) skip *= 2;
       else if ( fabs(ke_old - ke_sum)/ke_old  > 0.1 && skip > 1) skip /= 2;
     }
     MPI_Bcast(&skip, 1, MPI_INT, ROOT, MPI_COMM_WORLD);
  #else
     if ( fabs(ke_old - ke_sum)/ke_old < 0.0001) skip *= 2;
     else if ( fabs(ke_old - ke_sum)/ke_old  > 0.1 && skip > 1) skip /= 2;
  #endif
     ke_old = ke_sum;
     sv_sum = 0.0;
     ke_sum = 0.0;
     vz_sum = 0.0;
     for (ii=0; ii<np[0]; ii++)
     {
        energy = v_r[0][ii]*v_r[0][ii] + v_theta[0][ii]*v_theta[0][ii] + v_z[0][ii]*v_z[0][ii];
        ke_sum += energy;
        vz_sum += v_z[0][ii];
        vel = sqrt(energy);
        energy *= Escale[0];
	sv_sum += (*sigmaTptr)(energy)*vel;
     }
#ifdef MPI_1D
     sum[0] = sv_sum;
     sum[1] = vz_sum;
     MPI_Allreduce(sum, sum_tot, 2, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
     I_z = Iz_norm*sum_tot[1];
#else
     I_z =  Iz_norm*vz_sum;
#endif
  }
	    
#ifdef MPI_1D
  MPI_Allreduce(np, np_tot, NSMAX, MPI_INT, MPI_SUM, MPI_COMM_WORLD);
  Ne = (float) np_tot[0];
  mu_e = mu_norm*Ne/sum_tot[0];
  E_z = source()/(mu_e*Ne);
  /*if (my_rank == 0) printf("time=%e, E_z=%f, Iz=%f, Ne=%f, mu_e=%f\n", atime, E_z/(ptopn*dr), I_z, Ne, mu_e);*/
#else
  Ne = (float) np[0];
  mu_e = mu_norm*Ne/sv_sum;
  E_z = source()/(mu_e*Ne);
 /* printf("time=%e, E_z=%f, Iz=%f, Ne=%f, mu_e=%f\n", atime, E_z/(ptopn*dr), I_z, Ne, mu_e);*/
#endif
  yp[0] = rhon*rho[0]*vol[0];
  for (ii=1; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>=0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];

} /* end of circuit5 */
/***************************************************************/
/* 6: Case with no inner conductor and axial field supplied   */

void circuit6()
{
  int ii;
  float yp[nc];
  
  E_z = source();
  yp[0] = vol[0]*rhon*rho[0];
  for (ii=1; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>=0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];
}
/***************************************************************/
/* 7: Case with no inner conductor and axial power supplied   */

void circuit7()
{
  int ii;
  static int skip=1;
  float yp[nc];
  float energy, vz_sum;
#ifdef MPI_1D
  float vz_sum_tot;
#endif

  /* calculate I_z every skip timesteps */
  if (it[0]%skip == 0)
  {  
  #ifdef MPI_1D
     if (my_rank == ROOT) {
       if ( fabs(ke_old - ke_sum)/ke_old < 0.0001) skip *= 2;
       else if ( fabs(ke_old - ke_sum)/ke_old  > 0.1 && skip > 1) skip /= 2;
     }
     MPI_Bcast(&skip, 1, MPI_INT, ROOT, MPI_COMM_WORLD);
  #else
     if ( fabs(ke_old - ke_sum)/ke_old < 0.0001) skip *= 2;
     else if ( fabs(ke_old - ke_sum)/ke_old  > 0.1 && skip > 1) skip /= 2;
  #endif
     ke_old = ke_sum;
     ke_sum = 0.0;
     vz_sum = 0.0;
     for (ii=0; ii<np[0]; ii++)
     {
        energy = v_r[0][ii]*v_r[0][ii] + v_theta[0][ii]*v_theta[0][ii] + v_z[0][ii]*v_z[0][ii];
        ke_sum += energy;
        vz_sum += v_z[0][ii];
     }
#ifdef MPI_1D
     MPI_Allreduce(&vz_sum, &vz_sum_tot, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD);
     I_z = Iz_norm*vz_sum_tot;
#else
     I_z =  Iz_norm*vz_sum;
#endif
  }
  if (I_z) E_z = source()/I_z;
/*#ifdef MPI_1D
  if (my_rank == 0) printf("time=%e, E_z=%f, Iz=%f, Ne=%d\n", atime, E_z/(ptopn*dr), I_z, np[0]);
#else
   printf("time=%e, E_z=%f, Iz=%f, Ne=%d\n", atime, E_z/(ptopn*dr), I_z, np[0]);
#endif*/
  yp[0] = rhon*rho[0]*vol[0];
  for (ii=1; ii<nc; ii++) yp[ii]=xp[ii]*(rhon*rho[ii] + a_phi[ii]*yp[ii-1]);

  phi[nc] = 0.0;
  for (ii=nc-1; ii>=0; ii--) phi[ii] = yp[ii] + c_phi[ii]*xp[ii]*phi[ii+1];

}