start.c

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

  {
     sp_n[i]= (float  *)malloc(ng*sizeof(float));
     sp_n_k[i]= (float  *)malloc(ng*sizeof(float));
     sp_n_0[i]= (float  *)malloc(ng*sizeof(float));
     sp_n_mcc[i]= (float  *)malloc(ng*sizeof(float));
      
     for(j=0; j<ng; j++)
       sp_n[i][j] = sp_n_0[i][j] = sp_n_k[i][j] = sp_n_mcc[i][j] = 0.0;
  }
  
  rho= (float  *)malloc(ng*sizeof(float));
  rho_unfilt= (float  *)malloc(ng*sizeof(float));
  phi= (float  *)malloc(ng*sizeof(float));
  efield  = (float  *)malloc(ng*sizeof(float));
  phi2 = (float  *)malloc(ng*sizeof(float));
  efield2  = (float  *)malloc(ng*sizeof(float));

  for (j=0; j<ng; j++) rho_unfilt[j] = 0.0;
  
  /**** Allocate Dianostic arrays  *******/
  if(theRunWithXFlag)
  {
    jdote = (float **)malloc(nsp*sizeof(float *));
    jzdotez = (float **)malloc(nsp*sizeof(float *));
    j_r = (float **)malloc(nsp*sizeof(float *));
    jr_unfilt = (float **)malloc(nsp*sizeof(float *));
    jr_prof = (float **)malloc(nsp*sizeof(float *));
    j_z = (float **)malloc(nsp*sizeof(float *));
    jz_prof = (float **)malloc(nsp*sizeof(float *));
    mu_prof = (float **)malloc(nsp*sizeof(float *));
    ke_r = (float **)malloc(nsp*sizeof(float *));
    ker_r = (float **)malloc(nsp*sizeof(float *));
    ket_r = (float **)malloc(nsp*sizeof(float *));
    kez_r = (float **)malloc(nsp*sizeof(float *));
    ke_avg = (float **)malloc(nsp*sizeof(float *));
    ker_avg = (float **)malloc(nsp*sizeof(float *));
    ket_avg = (float **)malloc(nsp*sizeof(float *));
    kez_avg = (float **)malloc(nsp*sizeof(float *));
    for(isp=0; isp<nsp; isp++) 
    {    
       j_r[isp] = (float *)malloc(ng*sizeof(float));
       jr_unfilt[isp] = (float *)malloc(ng*sizeof(float));
       jr_prof[isp] = (float *)malloc(ng*sizeof(float));
       j_z[isp] = (float *)malloc(ng*sizeof(float));
       jz_prof[isp] = (float *)malloc(ng*sizeof(float));
       mu_prof[isp] = (float *)malloc(ng*sizeof(float));
       ke_r[isp] = (float *)malloc(ng*sizeof(float));
       ker_r[isp] = (float *)malloc(ng*sizeof(float));
       ket_r[isp] = (float *)malloc(ng*sizeof(float));
       kez_r[isp] = (float *)malloc(ng*sizeof(float));
       ke_avg[isp] = (float *)malloc(ng*sizeof(float));
       ker_avg[isp] = (float *)malloc(ng*sizeof(float));
       ket_avg[isp] = (float *)malloc(ng*sizeof(float));
       kez_avg[isp] = (float *)malloc(ng*sizeof(float));
       jdote[isp] = (float *)malloc(ng*sizeof(float));
       jzdotez[isp] = (float *)malloc(ng*sizeof(float));
       for(j=0; j<ng; j++) jr_unfilt[isp][j] = jdote[isp][j] = jzdotez[isp][j]
          = jr_prof[isp][j] = jz_prof[isp][j] = mu_prof[isp][j]
	  = ke_avg[isp][j] = ker_avg[isp][j] = ket_avg[isp][j] = kez_avg[isp][j]
	  = ke_r[isp][j] = ker_r[isp][j] = ket_r[isp][j] = kez_r[isp][j] = 0.0;
    }

  /* Allocate space for time history diagnostics */
  
    t_array= (float  *)malloc(HISTMAX*sizeof(float));
    ke_tot= (float  *)malloc(HISTMAX*sizeof(float));
    te_hist= (float  *)malloc(HISTMAX*sizeof(float));
    ese_hist= (float  *)malloc(HISTMAX*sizeof(float));
    sigma_hist= (float  *)malloc(HISTMAX*sizeof(float));
    Ez_hist= (float  *)malloc(HISTMAX*sizeof(float));
    mu_hist= (float  *)malloc(HISTMAX*sizeof(float));
    mu_hist2 = (float  *)malloc(HISTMAX*sizeof(float));
    mu_hist3 = (float  *)malloc(HISTMAX*sizeof(float));
    Iz_hist= (float  *)malloc(HISTMAX*sizeof(float));
    Iz_hist2= (float  *)malloc(HISTMAX*sizeof(float));
    nu_iz_hist= (float  *)malloc(HISTMAX*sizeof(float));
    nu_ex_hist= (float  *)malloc(HISTMAX*sizeof(float));
    nu_el_hist= (float  *)malloc(HISTMAX*sizeof(float));
  
    com_cur_hist= (float  *)malloc(HISTMAX*sizeof(float));
    com_phi_hist= (float  **)malloc(2*sizeof(float  *));
    com_phi_hist[0]= (float  *)malloc(HISTMAX*sizeof(float));
    com_phi_hist[1]= (float  *)malloc(HISTMAX*sizeof(float));
    sheath = (float **)malloc(2*sizeof(float *));
    sheath[0]=(float *)malloc(HISTMAX*sizeof(float));
    sheath[1]= (float *)malloc(HISTMAX*sizeof(float));
  
    com_power_hist = (float  *)malloc(HISTMAX*sizeof(float));
    avg_pow_hist = (float  *)malloc(HISTMAX*sizeof(float));
  
    np_hist= (float  **)malloc(nsp*sizeof(float  *));
    np_hist2= (float  **)malloc(nsp*sizeof(float  *));
    kes_hist= (float  **)malloc(nsp*sizeof(float  *));
    flux = (float  **)malloc(nsp*sizeof(float  *));

    C_hist = (float **)malloc(3*sizeof(float *));		/* LEEHJ */
    C_hist[0] = (float *)malloc(HISTMAX*sizeof(float));
    C_hist[1] = (float *)malloc(HISTMAX*sizeof(float));

    for(i=0; i<nsp; i++)
    {
      np_hist[i]= (float  *)malloc(HISTMAX*sizeof(float));
      np_hist2[i]= (float  *)malloc(HISTMAX*sizeof(float));
      kes_hist[i]= (float  *)malloc(HISTMAX*sizeof(float));
      flux[i]= (float  *)malloc(HISTMAX*sizeof(float));
    }
  /*--------------------------------------------------*/
  /* Allocate space for frequency history diagnostics and
  time-average diagnositics */
  
    if(nfft)
    {

/**** Profile and time-averaged profiles of diagnostic arrays ****/
      density= (float **)malloc(nsp*sizeof(float *));
      den_sum= (float **)malloc(nsp*sizeof(float *));
      den_ave= (float **)malloc(nsp*sizeof(float *));
      for(isp=0; isp<nsp; isp++)
      {
         density[isp] = (float *)malloc(ng*sizeof(float));
         den_sum[isp] = (float *)malloc(ng*sizeof(float));
         den_ave[isp] = (float *)malloc(ng*sizeof(float));
         for(i=0; i< ng; i++) den_sum[isp][i]= den_ave[isp][i]= 0.0;
      }
      phi_sum= (float *)malloc(ng*sizeof(float));
      phi_ave= (float *)malloc(ng*sizeof(float));
      for(i=0; i< ng; i++) phi_sum[i]= phi_ave[i]= 0.0;

      e_sum= (float *)malloc(ng*sizeof(float));
      e_ave= (float *)malloc(ng*sizeof(float));
      for(i=0; i< ng; i++) e_ave[i]= e_sum[i]= 0.0;

      rho_ave= (float *)malloc(ng*sizeof(float));
      rho_sum= (float *)malloc(ng*sizeof(float));
      for(i=0; i< ng; i++) rho_sum[i]= 0.0;

/**** Frequency history arrays ****/
      freq_hi= nfft/2;
      df= 1/(nfft*dt);

      f_array= (float  *)malloc(nfft*sizeof(float));
      for(i=0; i< nfft; i++) f_array[i]= i*df;
      cur_hist= (float  *)malloc(nfft*sizeof(float));
      phi_hist= (float  **)malloc(2*sizeof(float  *));
      phi_hist[0]= (float  *)malloc(nfft*sizeof(float));
      phi_hist[1]= (float  *)malloc(nfft*sizeof(float));
      pow_hist= (float  *)malloc(nfft*sizeof(float));
    
      Local_t_array = (float  *)malloc(nfft*sizeof(float));
    
/**** FFT arrays ****/
      cur_fft= (float  *)malloc(nfft*sizeof(float));
      phi_fft= (float  *)malloc(nfft*sizeof(float));
      mphi_fft=(float  *)malloc(nfft*sizeof(float));
      pow_fft= (float  *)malloc(nfft*sizeof(float));
      z_fft= (float  *)malloc(nfft*sizeof(float));
    }
  } /*end of if (theRunWithXFlag) loop */
  
  /*------------------------------------------*/
  /* Set up each "species" parameters         */
  vrscale=0.0;
  for (isp=0; isp<nsp; ++isp)
  {
    species(initnp,isp);
    
    /* Variables for converting from physical to code units */
    q_ratio[isp] = charge[isp]/ECHARGE;
    qm[isp]= charge[isp]/mass[isp];
    wc[isp] = qm[isp]*Bz*dt;
    Escale[isp] = 0.5*vscale*vscale*mass[isp]/ECHARGE;
    intoi[isp] = charge[isp]*nc2p/dt;      /* Normalisation for current */
    a_scale[isp] = qm[isp]*me_e*sp_k[isp];

    /* Renormalize all velocities */
    v0[isp][0]/= vscale;
    v0[isp][1]/= vscale;
    vc[isp][0]/= vscale;
    vc[isp][1]/= vscale;
    vt[isp][0]/= vscale;
    vt[isp][1]/= vscale;
    v0_th[isp] /= vscale;
    vt_th[isp] /= vscale;
    v0_z[isp] /= vscale;
    vt_z[isp] /= vscale;
    
    emin[isp] /= Escale[isp];
    de[isp]   /= Escale[isp];
    de_mid[isp]   /= Escale[isp];
    emin_mid[isp] /= Escale[isp];
   
    rs_mid[isp] /= dr;
    rf_mid[isp] /= dr;
  }
  
  if (RT_flag) input_for_rad();                          /* LEEHJ ; RAD-TR */
  else nc_RT=nc;
  ng_array=(float *)malloc(nc_RT*sizeof(float));
  for (j=0;j<nc_RT;j++) ng_array[j]=1.0;

  if (nc%nc_RT==0) nc_ratio=nc/nc_RT;
  else {
    puts("nc should be a multiple of nc_RT");
    exit(-1);
  }
  fclose(InputDeck);
  if (RT_flag) start_rad();

/* Normalise parameters which depend on nc2p (can change in species
sub-routine */
  rhoback /= dntod*ECHARGE;
  rhon = ptopn*(ECHARGE/EPS0)*nc2p/(PI*height);

  if(!vrscale) vrscale = 1e4;
  vrscale *= 5;

  /* Allocate particle arrays */
  
  r  = (float **)malloc(nsp*sizeof(float *));
  v_r = (float **)malloc(nsp*sizeof(float *));
  v_theta= (float **)malloc(nsp*sizeof(float *));
  v_z = (float **)malloc(nsp*sizeof(float *));
  angular_mtm = (float **)malloc(nsp*sizeof(float *));
  
  for(isp=0; isp< nsp; isp++)
  {
    r[isp]  = (float *)malloc(maxnp[isp]*sizeof(float));
    v_r[isp] = (float *)malloc(maxnp[isp]*sizeof(float));
    v_theta[isp]= (float *)malloc(maxnp[isp]*sizeof(float));
    v_z[isp] = (float *)malloc(maxnp[isp]*sizeof(float));
    angular_mtm[isp] = (float *)malloc(maxnp[isp]*sizeof(float));
    for(j=0; j<maxnp[isp]; j++)
      r[isp][j]=v_r[isp][j]=v_theta[isp][j]=v_z[isp][j]=angular_mtm[isp][j]=0.0;
  }
  
  if (nsp && (!r[nsp-1] || !v_r[nsp-1] || !v_theta[nsp-1] || !v_z[nsp-1]))
  {
  #ifdef MPI_1D
    if(my_rank == ROOT)
      puts("START: r, v malloc failed");
  #else
    puts("START: r, v malloc failed");
  #endif
    exit(1);
  }
   
  /*************************************************/
  /* Allocate space for the dist. functions arrays */
  
  if (theRunWithXFlag)
  {

/* Allocate arrays for velocity distribution diagnostic */
     if (vel_dist_accum)
     {
       vr_dist = (float ***)malloc(nsp*sizeof(float **));
       vth_dist = (float ***)malloc(nsp*sizeof(float **));
       vz_dist = (float ***)malloc(nsp*sizeof(float **));
       for(isp=0; isp<nsp; isp++)
       {
          vr_dist[isp] = (float **)malloc(ng*sizeof(float*));
          vth_dist[isp] = (float **)malloc(ng*sizeof(float*));
          vz_dist[isp] = (float **)malloc(ng*sizeof(float*));
          for (i=0; i<ng; i++)
          {
             vr_dist[isp][i] = (float *)calloc(nvrbin[isp],sizeof(float));
             vth_dist[isp][i] = (float *)calloc(nvtbin[isp],sizeof(float));
             vz_dist[isp][i] = (float *)calloc(nvzbin[isp],sizeof(float));
          }
       }
       vr_array = (float **)malloc(nsp*sizeof(float *));
       vth_array = (float **)malloc(nsp*sizeof(float *));
       vz_array = (float **)malloc(nsp*sizeof(float *));
       for(isp=0; isp<nsp; isp++)
       {
          vr_array[isp] = (float *)malloc(nvrbin[isp]*sizeof(float));
          vth_array[isp] = (float *)malloc(nvtbin[isp]*sizeof(float));
          vz_array[isp] = (float *)malloc(nvzbin[isp]*sizeof(float));
       }
     }

/* energy distribution functions */
    fe      = (float **)malloc(nsp*sizeof(float *));
    e_array = (float **)malloc(nsp*sizeof(float *));
    ftheta  = (float **)malloc(nsp*sizeof(float *));
    th_array= (float **)malloc(nsp*sizeof(float *));
    fe_mid = (float **)malloc(nsp*sizeof(float *));
    fe_mid_show = (float **)malloc(nsp*sizeof(float *));
    emid_array= (float **)malloc(nsp*sizeof(float *));
  
    for(isp=0; isp<nsp; isp++)
    {
      fe[isp]     = (float *)malloc(nbin[isp]*sizeof(float));
      ftheta[isp]  = (float *)malloc(nbin[isp]*sizeof(float));
      e_array[isp]= (float *)malloc(nbin[isp]*sizeof(float));
      th_array[isp]= (float *)malloc(nbin[isp]*sizeof(float));
      fe_mid[isp]  = (float *)malloc(nbin_mid[isp]*sizeof(float));
      fe_mid_show[isp]  = (float *)malloc(nbin_mid[isp]*sizeof(float));
      emid_array[isp]= (float *)malloc(nbin_mid[isp]*sizeof(float));
    }
  
    for(isp=0; isp<nsp; isp++)
    {
      for(i=0; i<nbin[isp]; i++)
      {
        fe[isp][i]      = 0.0;
        ftheta[isp][i]  = 0.0;
        e_array[isp][i] = emin[isp] +i*de[isp];
        th_array[isp][i]= i*dtheta[isp];
      }
      for (i=0; i<nbin_mid[isp]; i++)
      {
          emid_array[isp][i]= emin_mid[isp] +(i+.5)*de_mid[isp];
          fe_mid[isp][i] = fe_mid_show[isp][i] = 0.0;
      }
    }   
  }

#ifdef MPI_1D

  /* MPI parameters */
  counts = (int *)malloc(num_procs*sizeof(int));
  displs = (int *)malloc(num_procs*sizeof(int));

  /* MPI total arrays */
  sp_n_tot = (float **)malloc(nsp*sizeof(float *));
  for(isp=0; isp<nsp; isp++) {
      iwall_tot[isp] = 0;
      sp_n_tot[isp] = (float *)malloc(ng*sizeof(float));
      for(j=0; j<ng; j++)
         sp_n_tot[isp][j] = 0.0;