initwin.c

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

           sprintf(buffer, "V_th dist species %d", isp+1);
           XGSet3D("linlinlin", "R", "Vt", buffer, -100.0, 40.0, "closed",
                700, 175, dr, 1, 1.0, True, True, True, 0, 1.0, 0, 1.0, 0, 1.0);
           XGSurf(r_grid, vth_array[isp], vth_dist[isp], &ng, &nvtbin[isp], 4);
      }
      if (nvzbin[isp])
      {  
          sprintf(buffer, "V_z dist species %d", isp+1);
          XGSet3D("linlinlin", "R", "Vz", buffer, -100.0, 40.0, "closed",
                700, 175, dr, 1, 1.0, True, True, True, 0, 1.0, 0, 1.0, 0, 1.0);
          XGSurf(r_grid, vz_array[isp], vz_dist[isp], &ng, &nvzbin[isp], 4);
      }
    }
  } /* end if vel_dist_accum */

  /*********************************************/

  XGSet2D("linlin", "Theta", "f(theta)", "closed", 600, 350, 180/PI, 1.0,
          False, True, 0.0, 90.0, 0.0, 0.0);

  for (isp=nsp-1; isp>=0; isp--) XGCurve(th_array[isp], ftheta[isp],
         &nbin[isp], isp);

  /*********************************************/

  XGSet2D("linlog", "Time", "Ave Kinetic Energy(t) (eV)", "closed", 500, 200,
          1.0, 1.0, True, True, 0.0, 0.0, 0.0, 0.0);

  for (isp=nsp-1; isp>=0; isp--) XGCurve(t_array, kes_hist[isp], &hist_hi, isp);

  /*********************************************/

  XGSet2D("linlin", "Time", "Field Energy(t) (eV)", "closed", 650, 400, 1.0,
         1.0, True, True, 0.0, 0.0, 0.0, 0.0);

  XGCurve(t_array, ese_hist, &hist_hi, 4);
  XGCurve(t_array, ke_tot, &hist_hi, 5);
  XGCurve(t_array, te_hist, &hist_hi, 6);

  /*********************************************/

  XGSet2D("linlin", "Time", "Sheath: Live", "closed", 650, 400, 1.0, dr,
                True, False,0.0, 0.0, 0.0, (r1-r0)*dr);

  XGCurve(t_array, sheath[0], &hist_hi, 1);
/*  XGCurve(t_array, esheath[0], &hist_hi, 3); */

  /*********************************************/

  XGSet2D("linlin", "Time", "Sheath: Ground", "closed", 650, 400, 1.0, dr,
                True, False, 0.0, 0.0, 0.0, (r1-r0)*dr);

  XGCurve(t_array, sheath[1], &hist_hi, 1);
/*  XGCurve(t_array, esheath[1], &hist_hi, 3); */

  /*********************************************/
/* Collision profiles */
  if(ecollisional || icollisional)
  {
     for(ii=0; ii<ncolls; ii++)
     {
        XGSet2D("linlin", "R", rate_title[ii], "closed", 100, 500, 
          1.0, dntod/dt, False, True, r0*dr, r1*dr, 0.0, 1.0);

        XGCurve(r_array, rate_show[ii], &ng, 3);
     }
/*
     if (RT_flag) {
        XGSet2D("linlin", "R", "MCC Meta", "closed", 100, 500,
                Dr, dntod/dt, False, True, r0*dr, r1*dr, 0.0, 1.0);
        XGCurve(r_grid2, ex_rate_show[1], &nc_RT, 3);

        XGSet2D("linlin", "R", "MCC Rad", "closed", 100, 500,
                Dr, dntod/dt, False, True, r0*dr, r1*dr, 0.0, 1.0);
        XGCurve(r_grid2, ex_rate_show[0], &nc_RT, 3);
     }
*/
     XGSet2D("linlin", "R","Elastic Heat Loss", "closed",
              100, 500, 1.0, 1.0, False, True, r0*dr, r1*dr, 0.0, 1.0);

     XGCurve(r_array, Pel_show, &ng, 3);

     XGSet2D("linlin", "R","Inelastic Heat Loss", "closed",
              100, 500, 1.0, 1.0, False, True, r0*dr, r1*dr, 0.0, 1.0);

     XGCurve(r_array, Pie_show, &ng, 3);

     XGSet2D("linlin", "R","Charge Exchange Heat Loss", "closed",
              100, 400, 1.0, 1.0, False, True, r0*dr, r1*dr, 0.0, 1.0);

     XGCurve(r_array, Pcx_show, &ng, 3);
  }
  /*********************************************/

  if(nfft)
  {
    XGSet2D("linlin", "X", "Time Ave E field(r)", "close", 10, 500, 1.0, 
        1.0/(ptopn*dr), False, True, r0*dr,r1*dr, 0.0, 0.0);
   
    XGCurve(r_array, e_ave, &ng, 3);

    /*********************************************/

    XGSet2D("linlin", "X", "Time Ave Potential(r)", "close", 400, 300, 1.0, 
          1.0/ptopn, False, True, r0*dr,r1*dr, 0.0, 0.0);
   
    XGCurve(r_array, phi_ave, &ng, 4);

    /*********************************************/

    XGSet2D("linlin", "X", "Time Ave Density n(r)", "close", 600, 500, 1.0, 
	     1.0, False, True, r0*dr,r1*dr, 0.0, 0.0);
   
    for (isp=nsp-1; isp>=0; isp--) XGCurve(r_array, den_ave[isp], &ng, isp);

    /*********************************************/

   XGSet2D("linlin", "X", "Time Ave Charge n(r)", "close", 400, 300, 1.0, 1.0,
            False, True, r0*dr,r1*dr, 0.0, 0.0);
  
    XGCurve(r_array, rho_ave, &ng, 4);


    /*********************************************/

    XGSet2D("linlin", "Time", "Local Mid Potential(t)", "close", 300, 550,
            1.0, 1.0/ptopn, True, True, 0.0, 0.0, 0.0, 0.0);
   
    XGCurve(Local_t_array, phi_hist[1], &thist_hi, 3);

    /*********************************************/

    XGSet2D("linlin", "Time", "Local LHS Potential(t)", "close", 300, 550,
            1.0, 1.0/ptopn, True, True, 0.0, 0.0, 0.0, 0.0);

    XGCurve(Local_t_array, phi_hist[0], &thist_hi, 3);

    /*********************************************/

    XGSet2D("linlin", "Time", "Local Current(t)", "close", 300, 550, 1.0, 
    	    1.0, True, True, 0.0, 0.0, 0.0, 0.0);

    XGCurve(Local_t_array, cur_hist, &thist_hi, 3);

    /*********************************************/

    XGSet2D("linlog", "Freq", "Mag of Mid V(f)", "close", 100, 250, 1.0, 1.0, 
           True, True, 0.0, 0.0, 0.0, 0.0);

    XGCurve(f_array, mphi_fft, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlin", "Freq", "Phase of Mid V(f)", "close", 400, 250, 1.0,
             1.0, True, False, 1.0, 1.0, -180.0, 180.0);

    XGCurve(f_array, mphi_fft+freq_hi, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlog", "Freq", "Mag of LHS V(f)", "close", 100, 250, 1.0, 1.0, 
           True, True, 0.0, 0.0, 0.0, 0.0);

    XGCurve(f_array, phi_fft, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlin", "Freq", "Phase of LHS V(f)", "close", 400, 250, 1.0, 
            1.0, True, False, 1.0, 1.0, -180.0, 180.0);

    XGCurve(f_array, phi_fft+freq_hi, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlog", "Freq", "Mag of I(f)", "close", 100, 250, 1.0, 1.0,
        True, True, 0.0, 0.0, 0.0, 0.0);

    XGCurve(f_array, cur_fft, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlin", "Freq", "Phase of I(f)", "close", 400, 250, 1.0, 1.0,
            True, False, 1.0, 1.0, -180.0, 180.0);

    XGCurve(f_array, cur_fft+freq_hi, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlog", "Freq", "Mag of P(f)", "close", 100, 250, 1.0, 1.0,
        True, True, 0.0, 0.0, 0.0, 0.0);

    XGCurve(f_array, pow_fft, &freq_hi, 4);

    /*********************************************/

    XGSet2D("linlin", "Freq", "Phase of P(f)", "close", 400, 250, 1.0, 1.0,
            True, False, 1.0, 1.0, -180.0, 180.0);

    XGCurve(f_array, pow_fft+freq_hi, &freq_hi, 4);

    /*********************************************/
     XGSet2D("linlin", "Freq", "Phase btwn V & I", "close", 400, 250, 1.0, 1.0,
                    True, False, 1.0, 1.0, -180.0, 180.0);

      XGCurve(f_array, z_fft+freq_hi, &freq_hi, 4);

    /*********************************************/
  }
}
char Revision[]={'2','.','1'};
char RevisionR[]={'2','.','R'};

/***********************************************
** Dumps data to a file so that simulation can
** be restarted from the time of the dump
************************************************/
void Dump(char *filename)
{
  int i, j;
  float ftemp, length, temp;
  char dropfile[50];
  FILE *DMPFile, *DropFile;
 
#ifdef MPI_1D
  int isp, source, tag, sum; 
  MPI_Status status;
  
  if (my_rank == ROOT) { /* Only ROOT proc. dumps to the file */
    if ((DMPFile = fopen(filename, "w")) == NULL) {
      puts("Dump: open failed");
      return;
    }
    
    printf("Start Dump\n");
  }

  for (isp=0; isp < nsp; isp++){/* Gather data to dump into dumpfile */
      tag = 50 + isp;
      if(my_rank == ROOT)
      {
         counts[0] = np[isp];
         displs[0] = 0;
         sum = counts[0];
         for (source = 1; source <num_procs; source++){
           MPI_Recv(&counts[source], 1, MPI_INT, source,tag,
              MPI_COMM_WORLD, &status);
           displs[source] = displs[source - 1] + counts[source-1];
           sum = sum + counts[source];
         }
         np_tot[isp] = sum;
      }
      else{
        MPI_Send(&np[isp], 1, MPI_INT, 0, tag, MPI_COMM_WORLD);
      }
      
      MPI_Gatherv(r[isp], np[isp], MPI_FLOAT, r_tot[isp], counts, 
         displs, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
      MPI_Gatherv(v_r[isp], np[isp], MPI_FLOAT, vr_tot[isp], counts, 
         displs, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
      MPI_Gatherv(v_theta[isp], np[isp], MPI_FLOAT, vtheta_tot[isp], counts, 
         displs, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
      MPI_Gatherv(v_z[isp], np[isp], MPI_FLOAT, vz_tot[isp], counts, 
         displs, MPI_FLOAT, ROOT, MPI_COMM_WORLD);

     /* MPI_Barrier(MPI_COMM_WORLD);*/
  }

  if (my_rank == ROOT) { /* Only ROOT proc. dumps to the file */
 
    printf("np_tot[0] = %d\n", np_tot[0]);
    printf("proc 0:  np[0] = %d\n", np[0]);
    if(nsp > 1)
      {
	printf("np_tot[1] = %d\n", np_tot[1]);
	printf("proc 0:  np[1] = %d\n", np[1]);
      }

    printf("sp_n_tot[0][ng/2] = %g\n", sp_n_tot[0][ng/2]);
    if(nsp > 1)
      {
	printf("sp_n_tot[1][ng/2] = %g\n", sp_n_tot[1][ng/2]);
      }
    printf("proc 0: sp_n[0][ng/2] = %g\n", sp_n[0][ng/2]/vol[ng/2]);
    if(nsp > 1)
      {
	printf("proc 0: sp_n[1][ng/2] = %g\n", sp_n[1][ng/2]/vol[ng/2]);
      }

    if (RT_flag) XGWrite(RevisionR, 1, 4, DMPFile, "char");
    else XGWrite(Revision, 1, 4, DMPFile, "char");
    
    ftemp = atime;
    XGWrite(&ftemp, 4, 1, DMPFile, "float");
    XGWrite(&sigma, 4, 1, DMPFile, "float");
    XGWrite(&q_0, 4, 1, DMPFile, "float");
    XGWrite(&q_1, 4, 1, DMPFile, "float");
    XGWrite(&q_2, 4, 1, DMPFile, "float");
    XGWrite(&q_3, 4, 1, DMPFile, "float");
    XGWrite(&nsp, 4, 1, DMPFile, "int");
    XGWrite(iwall_tot, 4, nsp, DMPFile, "float");
    XGWrite(np_tot, 4, nsp, DMPFile, "int");
    XGWrite(&nc2p, 4, 1, DMPFile, "float");
    
    length = r1 - r0;
    for (isp=0; isp<nsp; isp++) {
      for (i=0; i<np_tot[isp]; i++) {
        ftemp = (r_tot[isp][i] -r0)/length;
        XGWrite(&ftemp, 4, 1, DMPFile, "float");
        ftemp = vscale*vr_tot[isp][i];
        XGWrite(&ftemp, 4, 1, DMPFile, "float");
        ftemp = vscale*vtheta_tot[isp][i];
        XGWrite(&ftemp, 4, 1, DMPFile, "float");
        ftemp = vscale*vz_tot[isp][i];
        XGWrite(&ftemp, 4, 1, DMPFile, "float");
      }
    }
    /*------ Save data for Radiation Transport --------*/

    if (RT_flag) {
      for (i=0;i<nc_RT;i++) {
        ftemp=sp_ex[i];
        XGWrite(&ftemp, 4, 1, DMPFile, "float");
        ftemp=sp_exm[i];
        XGWrite(&ftemp, 4, 1, DMPFile, "float");
      }
      for (i=0;i<2;i++) {
        for (j=0;j<nc_RT;j++) {
          ftemp=ex_rate_show[i][j];
          XGWrite(&ftemp, 4, 1, DMPFile, "float");
        }
      }
      ftemp=(float) it_rate;
      XGWrite(&ftemp, 4, 1, DMPFile, "float");
    }

    fclose(DMPFile);
    /***  Then write the number of particles and Ez in ascii  *****/
    strcpy(dropfile, filename);
    strcat(dropfile, ".drp");

    DropFile = fopen(dropfile, "a");
    fprintf(DropFile, "%e ", atime);
    for(i=0; i<nsp; i++) fprintf(DropFile, "%d ", np_tot[i]);
    fprintf(DropFile, "%f ", E_z/(ptopn*dr));
    fprintf(DropFile, "%f ", I_z);
/* 
    if (RT_flag) {
      length=ftemp=0;
      for (i=0;i<nc_RT;i++) {
        temp=2.0*i+1;
        length+=sp_ex[i]*temp;
        ftemp+=sp_exm[i]*temp;
      }
      temp=(float)(nc_RT); temp*=temp;
      length*=dntod/temp;