radtr.c

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

#include "pdc1.h"

void  calculate_Akm();
float Ki2(float), Ki1(float), zK1(float), IB1(float), omega(float),
      ei1m(float), GE0(float), J0(float);
float aK1I0(float, float), bK0I1(float, float), I1(float),
      phi_integral1(int, int, float), phi_integral2(int, int, float);

float dfreqx, **BB, *Diff, *Diffh, *r_grid2h;
float Kr=2e-13, Kcq=3e-21, Kiza=6.2e-16;

float I1_l[7]={ 0.5,        0.87890594, 0.51498869, 0.15084934, 
                0.02658733, 0.00301532, 0.00032411 },
      I1_u[9]={ 0.39894228,-0.03988024,-0.00362018, 0.00163801, -0.01031555, 
                0.02282967,-0.02895312, 0.01787654,-0.00420059 },
      K0_l[7]={-0.57721566, 0.42278420, 0.23069756, 0.03488590,
                0.00262698, 0.00010750, 0.00000740 },
      K0_u[7]={ 1.25331414,-0.07832358, 0.02189568,-0.01062446,
                0.00587872,-0.00251540, 0.00053208 },
      I0_l[7]={ 1.0,        3.5156229,  3.0899424,  1.2067492, 
                0.2659732,  0.0360768,  0.0045813 },
      I0_u[9]={ 0.39894228, 0.01328592, 0.00225319,-0.00157565, 0.00916281,
               -0.02057706, 0.02635537,-0.01647633, 0.00392377 },
      K1_l[7]={ 1.0,        0.15443144,-0.67278579,-0.18156897, 
               -0.01919402,-0.00110404,-0.00004686 },
      K1_u[7]={ 1.25331414, 0.23498619,-0.03655620, 0.01504268,
               -0.00780353, 0.00325614,-0.00068245 },
      K0uI1u[15], K0uI1l[19], I0lI1l[12], K0lI1l[12], 
      K1uI0u[15], K1uI0l[19], I1lI0l[12], K1lI0l[12];

 /************** Initialize Excited State Calculation ***************/
void start_rad()
{
  register int i,j;
  float dt_drdr, dt_tau, F;
  float epsi_k, sigma, sigmasq, temp_in_K;
  char filename[20], a_char[80];
  FILE *InputDeck;

  Dr = dr*nc_ratio;
  dt_drdr = dt*sp_k_ex/Dr/Dr;
  dt_tau  = dt*sp_k_ex*A_ki;
#ifdef MPI_1D
  if (my_rank == ROOT) printf("dt/tau=%g\n",dt_tau);
#else
  printf("dt/tau=%g\n",dt_tau);
#endif
  dfreqx  = (xmax-xmin)/xbin;
  xbin++;
  interval2=1;			    /* It is for history of t-dep. variables */
  I_factor*=dntod;
  Kr*=dt*dntod;
  Kcq*=dt;
  Kiza*=dt*dntod;

  r_grid2 = (float *)malloc(nc_RT*sizeof(float));
  r_grid2h= (float *)malloc((nc_RT+1)*sizeof(float));
  for(j=0; j<nc_RT; j++) r_grid2[j] = (float)j+0.5;
  for(j=0; j<(nc_RT+1); j++) r_grid2h[j]= (float)j;
  A     = (float **)malloc(nc_RT*sizeof(float *));
  AA    = (float **)malloc(nc_RT*sizeof(float *));
  BB    = (float **)malloc(3 *sizeof(float *));
  vol2  = (float *)malloc(nc_RT*sizeof(float));
  sp_ex = (float *)malloc(nc_RT*sizeof(float));
  sp_exm= (float *)malloc(nc_RT*sizeof(float));
  prod  = (float *)malloc(nc_RT*sizeof(float));
  prodm = (float *)malloc(nc_RT*sizeof(float));
  Diff  = (float *)malloc(nc_RT*sizeof(float));
  Diffh = (float *)malloc((nc_RT+1)*sizeof(float));
  error = (float *)malloc(nc_RT*sizeof(float));
  E_prod = (float **)malloc(4*sizeof(float *));
  ex_ave = (float *)malloc(nc_RT*sizeof(float));
  exm_ave = (float *)malloc(nc_RT*sizeof(float));
  eden_ave = (float *)malloc(nc_RT*sizeof(float));
  ex_mccrate = (float **)malloc(2*sizeof(float *));
  ex_ave_show = (float *)malloc(nc_RT*sizeof(float));
  ex_rate_show = (float **)malloc(2*sizeof(float *));
  exm_ave_show = (float *)malloc(nc_RT*sizeof(float));
  emission_coeff = (float *)malloc(nc_RT*sizeof(float));
  radiation_flux = (float *)malloc(nc_RT*sizeof(float));
  error_sum=(float *)malloc(HISTMAX*sizeof(float));
  t_array3 =(float *)malloc(HISTMAX*sizeof(float));
  tot_intensity=(float *)malloc(HISTMAX*sizeof(float));
  for (i=0;i<nc_RT;i++) error[i]=0;
  for (i=0;i<HISTMAX;i++) error_sum[i]=0;
  for (i=0;i<4;i++) E_prod[i]=(float *)malloc(nc_RT*sizeof(float));

/*  for (i=0;i<nc_RT;i++)
    ng_array[i]=(C1*(i+0.5)*(nc_RT-i-0.5) + N2)/temp;	*/

  I_ave_lhs = (float *)malloc(xbin*sizeof(float));
  I_ave_rhs = (float *)malloc(xbin*sizeof(float));
  I_ave_lhs_show = (float *)malloc(xbin*sizeof(float));
  I_ave_rhs_show = (float *)malloc(xbin*sizeof(float));
  x_array=(float *)malloc(xbin*sizeof(float));
  L_array=(float *)malloc(xbin*sizeof(float));
  intensity_lhs=(float *)malloc(xbin*sizeof(float));
  intensity_rhs=(float *)malloc(xbin*sizeof(float));

  t_array2=(float *)malloc(HISTMAX2*sizeof(float));
  st_nex=(float **)malloc(HISTMAX2*sizeof(float *));
  st_nexm=(float **)malloc(HISTMAX2*sizeof(float *));
  st_I_lhs=(float **)malloc(HISTMAX2*sizeof(float *));
  st_I_rhs=(float **)malloc(HISTMAX2*sizeof(float *));
  for(i=0; i< HISTMAX2; i++) {
    st_nex[i]=(float *)malloc(nc_RT*sizeof(float));
    st_nexm[i]=(float *)malloc(nc_RT*sizeof(float));
    st_I_lhs[i]=(float *)malloc(xbin*sizeof(float));
    st_I_rhs[i]=(float *)malloc(xbin*sizeof(float));
  }
  for (i=0;i<nc_RT;i++) {
    A [i] = (float *)malloc(nc_RT*sizeof(float));
    AA[i] = (float *)malloc(nc_RT*sizeof(float));
    sp_ex[i]=sp_exm[i]=prod[i]=prodm[i]=0;
    vol2[i]=(2*i+1)*nc_ratio*nc_ratio;
  }
  for (i=0;i<3;i++) BB[i] = (float *)malloc(nc_RT*sizeof(float));
  for (i=0;i<2;i++) {
    ex_mccrate[i] = (float *)malloc(nc_RT*sizeof(float));
    ex_rate_show[i] = (float *)malloc(nc_RT*sizeof(float));
  }
  switch (lineshape) {
    case LORENTZIAN :
         k_function=lorentz;
         F_factor=I_factor*delta_nu*0.5/Cx*k0;
         break;
    case DOPPLER :
         k_function=doppler;
         F_factor=I_factor*delta_nu/(2*sqrt(log(2)))/Cx*k0;
         break;
    case VOIGT :
         k_function=voigt;
         break;
  }
  for (i=0;i<xbin;i++) {
    x_array[i]=xmin+i*dfreqx;
    L_array[i]=k_function(x_array[i]);
  }
 /************ Diffusion Coefficient **************/
  epsi_k=124.0;				/* for Ar */
  sigma=3.418;				/* for Ar */
  sigmasq=sigma*sigma;
  temp_in_K=11609*T_gas;
  for (i=0;i<nc_RT;i++) {   
    /* For Ar gas */
    Diff[i]=1.997279e-4*sqrt(temp_in_K/39.948)*temp_in_K/
            (pressure*sigmasq*omega(temp_in_K/epsi_k));		
  }   
#ifdef MPI_1D
  if (my_rank == ROOT) 
    printf("pressure=%g Torr, diff_coeff.=%g m^2/s\n", pressure, Diff[3]);
#else
  printf("pressure=%g Torr, diff_coeff.=%g m^2/s\n", pressure, Diff[3]);
#endif
  for (i=1;i<nc_RT;i++) Diffh[i]=0.5*(Diff[i]+Diff[i-1]);
  Diffh[0]=0;
  Diffh[nc_RT]=Diffh[nc_RT-1];

 /************ Matrix element A_k,m ***************/
#ifdef MPI_1D
 if (my_rank == ROOT) {
  printf("The file name for A_km:"); scanf("%s",filename);
  if ((InputDeck = fopen(filename, "r"))) {
    puts("%%%%%%% Read A_km elements from the file %%%%%%%%");
    while(fscanf(InputDeck,"%f",&F)<1) fscanf(InputDeck,"%s",a_char);
    for (i=0;i<nc_RT;i++)
      for (j=0;j<nc_RT;j++)
        fscanf(InputDeck,"%f %f %f", &F, &F, &A[i][j]);
    fclose(InputDeck);
  }
 }
#else
  printf("The file name for A_km:"); scanf("%s",filename);
  if ((InputDeck = fopen(filename, "r"))) {
    puts("%%%%%%% Read A_km elements from the file %%%%%%%%");
    while(fscanf(InputDeck,"%f",&F)<1) fscanf(InputDeck,"%s",a_char);
    for (i=0;i<nc_RT;i++)
      for (j=0;j<nc_RT;j++)
        fscanf(InputDeck,"%f %f %f", &F, &F, &A[i][j]);
    fclose(InputDeck);
  }
#endif
  else {
#ifdef MPI_1D
   if (my_rank == ROOT) {
    printf("read A_km(): file open failed on %s.\n", filename);
    puts("%%%%%%% Calculate A_km elements %%%%%%%%");
   }
#else
    printf("read A_km(): file open failed on %s.\n", filename);
    puts("%%%%%%% Calculate A_km elements %%%%%%%%");
#endif
    calculate_Akm();
  }

 /******** Matrix element AA_k,m and BB_k,m *******/
  for (i=0;i<nc_RT;i++) {
    for (j=0;j<nc_RT;j++) AA[i][j] = dt_tau*A[i][j];
    AA[i][i]-= dt_tau;
    BB[1][i] = - dt_drdr*
               (Diffh[i]*r_grid2h[i]+Diffh[i+1]*r_grid2h[i+1])/r_grid2[i];
    if (i>0)     BB[0][i] = Diffh[i]*r_grid2h[i]/r_grid2[i]*dt_drdr;
    if (i<nc_RT) BB[2][i] = Diffh[i+1]*r_grid2h[i+1]/r_grid2[i]*dt_drdr;
  }
  BB[0][0]=0; BB[2][nc_RT-1]=0;

}
/*******************************************************************/

/*-------------- Evolution of Excited State --------------*/
void excited_evolution()
{
  register int i,j;
  static   int init_flag=1, count3=1, interval3=1;
  static float *temp, *tempm, *ei_pair, *half_ex, *half_exm, vel;
  float  Kr_Ne_Nm, Kcq_Ng, Kiza_ex, Kiza_mt;
  int    N_created, no;

  if (init_flag) {
    temp=(float *)malloc(nc_RT*sizeof(float));
    tempm=(float *)malloc(nc_RT*sizeof(float));
    half_ex=(float *)malloc(nc_RT*sizeof(float));
    half_exm=(float *)malloc(nc_RT*sizeof(float));
    ei_pair = (float *)malloc(nc_RT*sizeof(float));
    init_flag=0;
    vel=1.6241e6;
    hist_hi3=0;
    for (i=0;i<nc_RT;i++) ei_pair[i]=0;
  }

  for (i=0;i<nc_RT;i++) {
    temp[i]=sp_ex[i];
    tempm[i]=sp_exm[i];
  }

  /*  for (i=0; i< nsmoothing; i++) One_2_One(prod, nc_RT-1); 	*/
  /*-------------- Two step calculation ------------------*/
  for (i=0;i<nc_RT;i++) {
    Kr_Ne_Nm=0.5*Kr*eden_ave[i]*tempm[i];
    Kcq_Ng  =0.5*Kcq*gas_den;
    Kiza_ex =Kiza*temp[i]*temp[i];
    Kiza_mt =Kiza*tempm[i]*tempm[i];
    half_ex[i] =temp[i]+prod[i]*0.5;
    half_exm[i]=tempm[i]+prodm[i]*0.5;
    for (j=0;j<nc_RT;j++) half_ex[i]+=0.5*AA[i][j]*temp[j];
    half_ex[i]+=0.5*BB[1][i]*temp[i]; 
    half_exm[i]+=0.5*BB[1][i]*tempm[i];
    if (i>0) {
      half_ex[i]+=0.5*BB[0][i]*temp[i-1]; 
      half_exm[i]+=0.5*BB[0][i]*tempm[i-1];
    }
    if (i<nc_RT-1) {
      half_ex[i]+=0.5*BB[2][i]*temp[i+1]; 
      half_exm[i]+=0.5*BB[2][i]*tempm[i+1];
    }
    half_ex[i] += Kr_Ne_Nm - Kcq_Ng*temp[i] - Kiza_ex;
    half_exm[i]-= Kr_Ne_Nm + Kcq_Ng*tempm[i] + Kiza_mt;
  }
/*  for (i=0;i<nc_RT;i++) half_ex[i]=J0(2.405*(i+0.5)/nc_RT); */
  for (i=0;i<nc_RT;i++) {
    Kr_Ne_Nm=Kr*eden_ave[i]*half_exm[i];
    Kcq_Ng  =Kcq*gas_den;             
    Kiza_ex = Kiza*half_ex[i]*half_ex[i];
    Kiza_mt = Kiza*half_exm[i]*half_exm[i];
    E_prod[1][i]=prod[i]; 
    E_prod[2][i]=E_prod[3][i]=0;

    sp_ex[i]= temp[i]+prod[i];
    sp_exm[i]=tempm[i]+prodm[i]; 
    prod[i]=prodm[i]=eden_ave[i]=0;
    for (j=0;j<nc_RT;j++) sp_ex[i]+=AA[i][j]*half_ex[j];
    sp_ex[i] +=BB[1][i]*half_ex[i]; 
    sp_exm[i]+=BB[1][i]*half_exm[i];
    E_prod[2][i]+=BB[1][i]*half_ex[i]; 
    if (i>0) {
      sp_ex[i] +=BB[0][i]*half_ex[i-1];
      sp_exm[i]+=BB[0][i]*half_exm[i-1];
      E_prod[2][i]+=BB[0][i]*half_ex[i-1];
    }
    if (i<nc_RT-1) {
      sp_ex[i] +=BB[2][i]*half_ex[i+1];
      sp_exm[i]+=BB[2][i]*half_exm[i+1];
      E_prod[2][i]+=BB[2][i]*half_ex[i+1];
    }
/*    else {
      sp_ex[i] +=BB[2][i]*(2*half_ex[i]-half_ex[i-1]);
      sp_exm[i]+=BB[2][i]*(2*half_exm[i]-half_ex[i-1]);
      E_prod[2][i]+=BB[2][i]*(2*half_ex[i]-half_ex[i-1]);
    } */
    sp_ex[i] += Kr_Ne_Nm - Kcq_Ng*half_ex[i]  - 2*Kiza_ex;
    sp_exm[i]-= Kr_Ne_Nm + Kcq_Ng*half_exm[i] + 2*Kiza_mt;
    E_prod[3][i]+= Kr_Ne_Nm - Kcq_Ng*half_ex[i]  - 2*Kiza_ex;  
    ei_pair[i]+=(Kiza_ex+Kiza_mt)*vol2[i];
    E_prod[0][i]=E_prod[1][i]+E_prod[2][i]+E_prod[3][i];
  }

 /********** A* + A* -> e + A + A^+ ***********/
  for (i=0;i<nc_RT;i++) {
    N_created=0;
    while (ei_pair[i]>=nc2p) {
      N_created++;
      ei_pair[i]-=nc2p;
    }
    for (j=0;j<N_created;j++) { 
      no=np[0];
      r[0][no]=i+(j+1.0)/(N_created+1.0);
      r[0][no]*=nc_ratio;
      maxwellv(&v_r[0][no], &v_theta[0][no], &v_z[0][no], vel);
      no=np[1];
      r[1][no]=i+(j+1.0)/(N_created+1.0);
      r[1][no]*=nc_ratio;
      maxwellv(&v_r[1][no], &v_theta[1][no], &v_z[1][no], vg_therm);
      
   #ifdef MPI_1D
      if (my_rank == ROOT) {
        printf("electron: Vr=%g, Vth=%g, Vz=%g\n",
              v_r[0][np[0]],v_theta[0][np[0]],v_z[0][np[0]]);
        printf("ion: Vr=%g, Vth=%g, Vz=%g\n",
              v_r[1][np[1]],v_theta[1][np[1]],v_z[1][np[1]]);
      }
   #else
      printf("electron: Vr=%g, Vth=%g, Vz=%g\n",
              v_r[0][np[0]],v_theta[0][np[0]],v_z[0][np[0]]);
      printf("ion: Vr=%g, Vth=%g, Vz=%g\n",
              v_r[1][np[1]],v_theta[1][np[1]],v_z[1][np[1]]);
   #endif
   
      if (++np[0]>maxnp[0] || ++np[1]>maxnp[1]) {
      #ifdef MPI_1D
        if (my_rank == ROOT) 
	  puts("ADJUST(E-I PAIRS): too many particles. MUST EXIT!");
      #else
        puts("ADJUST(E-I PAIRS): too many particles. MUST EXIT!");
      #endif
          exit(-1);
      }
    }
  }  
 /************ Radiation Intensity ************/
  for (i=0;i<nc_RT;i++) {
    emission_coeff[i]=sp_ex[i];
    for (j=0;j<nc_RT;j++) emission_coeff[i]-=A[i][j]*sp_ex[j];
  }
  for (i=0;i<nc_RT;i++) {
    radiation_flux[i]=emission_coeff[i]*(i+0.25);
    for (j=0;j<i;j++) radiation_flux[i]+=emission_coeff[j]*(2*j+1.0);
    radiation_flux[i]/=(i+0.5);  	/* 2PI*Dr is multiplied in initwin.c */
  }

  /******** Variation of excited state ********/
  for (i=0;i<nc_RT;i++) {
    error[i]=E_prod[0][i];
    for (j=0;j<nc_RT;j++) error[i]+=AA[i][j]*sp_ex[j];
    if (sp_ex[i]) error[i]/=sp_ex[i];
    else error[i]=0;
  }

  /************ Time History ************/
  if (--count3==0) {
    if (hist_hi3>=HISTMAX) {
      for (i=1, j=4; i<HISTMAX/4; i++, j+=4) {
        error_sum[i]=error_sum[j];
        t_array3[i]=t_array3[j];
        tot_intensity[i]=tot_intensity[j];
      }
      hist_hi3=i;
      interval3*=4;
    }
    t_array3[hist_hi3]=atime;
    error_sum[hist_hi3]=0;
    for (i=0;i<nc_RT;i++) error_sum[hist_hi3]+=fabs(error[i]);
    tot_intensity[hist_hi3]=radiation_flux[nc_RT-1];
    hist_hi3++;

    count3 = interval3;
  }

}
/*******************************************************************/