initwin.c

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

      ftemp *=dntod/temp;
      for(i=0; i<nsp; i++) fprintf(DropFile, "%e ",np_tot[i]*nc2p/reactor_vol);
      fprintf(DropFile, "%e %e ", length, ftemp);
    }
*/
    fprintf(DropFile, "\n");
    fclose(DropFile);
    printf("Dump done\n");
  }
#else
  if ((DMPFile = fopen(filename, "w")) == NULL)
  {
    puts("Dump: open failed");
    return;
  }
  
  printf("Start Dump\n");
  printf("np[0] = %d\n", np[0]);
  printf("np[1] = %d\n", np[1]);

  printf("sp_n[0][ng/2] = %g\n", sp_n[0][ng/2]);
  printf("sp_n[1][ng/2] = %g\n", sp_n[1][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, 4, nsp, DMPFile, "float");
  XGWrite(np, 4, nsp, DMPFile, "int");
  XGWrite(&nc2p, 4, 1, DMPFile, "float");
  
  length = r1 - r0;
  for (j=0; j<nsp; j++)
  {
    for (i=0; i<np[j]; i++)
    {
      ftemp = (r[j][i] -r0)/length;
      XGWrite(&ftemp, 4, 1, DMPFile, "float");
      ftemp = v_r[j][i]*vscale;
      XGWrite(&ftemp, 4, 1, DMPFile, "float");
      ftemp = v_theta[j][i]*vscale;
      XGWrite(&ftemp, 4, 1, DMPFile, "float");
      ftemp = v_z[j][i]*vscale;
      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[i]);
  fprintf(DropFile, "%f ", E_z/(ptopn*dr));
  fprintf(DropFile, "%f ", I_z);
  fprintf(DropFile, "\n");
  fclose(DropFile);

  printf("Dump done\n");
#endif
}

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

void Restore(char *filename)
{
  char Rev[5];
  int i, j;
  float ftemp, length;
  FILE *DMPFile;

#ifdef MPI_1D

  float temp;
  int   isp, index, N, nnp, dest, tag, count;
  MPI_Status status;

  if (my_rank == ROOT){/* ROOT reads the dumpfile and broadcasts info. */

     if ((DMPFile = fopen(filename, "r+b")) == NULL) {
        puts("Dump: open failed");
        return;
     }

     XGRead(Rev, 1, 4, DMPFile, "char");
     for (i=0; i<4; i++)
       if (Rev[i]!=Revision[i] && Rev[i]!=RevisionR[i]) {
         puts("Incompatible dump file version");
         exit(1);
       }

     XGRead(&ftemp, 4, 1, DMPFile, "float");
     atime = ftemp;
     XGRead(&sigma, 4, 1, DMPFile, "float");
     XGRead(&q_0, 4, 1, DMPFile, "float");
     XGRead(&q_1, 4, 1, DMPFile, "float");
     XGRead(&q_2, 4, 1, DMPFile, "float");
     XGRead(&q_3, 4, 1, DMPFile, "float");
     XGRead(&nsp, 4, 1, DMPFile, "int");
     XGRead(iwall_tot, 4, nsp, DMPFile, "float");
     XGRead(np_tot, 4, nsp, DMPFile, "int");
     XGRead(&nc2p, 4, 1, DMPFile, "float");

     printf("Reading Restore:\n");
     printf("time = %g\n", atime);
     printf("nsp = %d\n", nsp);
     printf("np_tot[0] = %d\n", np_tot[0]);
     printf("np_tot[1] = %d\n", np_tot[1]);

  } /* end if (my_rank == ROOT) */

  /* Send parameter info to other procs */
  MPI_Bcast(&atime, 1, MPI_DOUBLE, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&sigma, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&q_0, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&q_1, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&q_2, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&q_3, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&iwall_tot, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&np_tot, nsp, MPI_INT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&nsp, 1, MPI_INT, ROOT, MPI_COMM_WORLD);
  MPI_Bcast(&nc2p, 1, MPI_FLOAT, ROOT, MPI_COMM_WORLD);

  oldsigma=sigma;   /* This is a cludge until I can put oldsigma into dump */
  length = r1 - r0;
     
  for (isp=0; isp<nsp; isp++) {
   if(my_rank == ROOT){
     for (i=0; i<np_tot[isp]; i++) {
        XGRead(&ftemp, 4, 1, DMPFile, "float");
        r_tot[isp][i] = ftemp*length + r0;
        XGRead(&ftemp, 4, 1, DMPFile, "float");
        vr_tot[isp][i]= ftemp/vscale;
        XGRead(&ftemp, 4, 1, DMPFile, "float");
        vtheta_tot[isp][i]= ftemp/vscale;
        XGRead(&ftemp, 4, 1, DMPFile, "float");
        vz_tot[isp][i]= ftemp/vscale;
     }

      /* Divvy up particles and send to various procs */
      nnp = np_tot[isp];
      np[isp] = floor(((float) np_tot[isp])/num_procs + 0.5);
      N = np[isp];
      printf("proc 1 to proc %d: isp = %d; np = %d\n",num_procs-1, isp, np[isp]);

      for(dest = 1; dest < num_procs; dest++){
        /* Randomly select N particles and move them to end of array */
        for(j=0; j< N; j++) {
           index= nnp*frand();
           nnp--;
           temp = r_tot[isp][nnp];
           r_tot[isp][nnp] = r_tot[isp][index];
           r_tot[isp][index] = temp;

           temp = vr_tot[isp][nnp];
           vr_tot[isp][nnp] = vr_tot[isp][index];
           vr_tot[isp][index] = temp;

           temp = vtheta_tot[isp][nnp];
           vtheta_tot[isp][nnp] = vtheta_tot[isp][index];
           vtheta_tot[isp][index] = temp;

           temp = vz_tot[isp][nnp];
           vz_tot[isp][nnp] = vz_tot[isp][index];
           vz_tot[isp][index] = temp;
        }

        for(j=0; j<N; j++) {
           r[isp][j] = r_tot[isp][nnp+j];
           v_r[isp][j] = vr_tot[isp][nnp+j];
           v_theta[isp][j] = vtheta_tot[isp][nnp+j];
           v_z[isp][j] = vz_tot[isp][nnp+j];
        }

        tag = 10 + isp;
        MPI_Send(&np[isp],1, MPI_INT,dest,tag,MPI_COMM_WORLD);
        tag = 100 + isp;
        count = np[isp];
        MPI_Send(r[isp], count, MPI_FLOAT,dest,tag,MPI_COMM_WORLD);
        tag = 200 + isp;
        MPI_Send(v_r[isp],count, MPI_FLOAT,dest,tag,MPI_COMM_WORLD);
        tag = 300 + isp;
        MPI_Send(v_theta[isp],count, MPI_FLOAT,dest,tag,MPI_COMM_WORLD);
        tag = 400 + isp;
        MPI_Send(v_z[isp],count, MPI_FLOAT,dest,tag,MPI_COMM_WORLD);

     }/* ends for(dest...) */

     /* Rest of particles assigned to ROOT proc. */
     for(j=0; j<nnp; j++) {
         r[isp][j] = r_tot[isp][j];
         v_r[isp][j] = vr_tot[isp][j];
         v_theta[isp][j] = vtheta_tot[isp][j];
         v_z[isp][j] = vz_tot[isp][j];
     }
     np[isp] = nnp; /* Note: nnp does not have to equal N */

     printf("proc 0: isp = %d; np = %d\n", isp, np[isp]);

   }/* ends if(my_rank == ROOT) */
   else{/* my_rank != ROOT */
     tag = 10 + isp;
     MPI_Recv(&np[isp],1, MPI_INT,ROOT,tag,MPI_COMM_WORLD,&status);
     tag = 100 + isp;
     count = np[isp];
     MPI_Recv(r[isp], count, MPI_FLOAT,ROOT,tag,MPI_COMM_WORLD,&status);
     tag = 200 + isp;
     MPI_Recv(v_r[isp],count, MPI_FLOAT,ROOT,tag,MPI_COMM_WORLD,&status);
     tag = 300 + isp;
     MPI_Recv(v_theta[isp],count, MPI_FLOAT,ROOT,tag,MPI_COMM_WORLD,&status);
     tag = 400 + isp;
     MPI_Recv(v_z[isp],count, MPI_FLOAT,ROOT,tag,MPI_COMM_WORLD,&status);
   }
 }/* ends for(isp...) */

 if (my_rank == ROOT) {
   /*------ Load data for Radiation Transport --------*/
     if (Rev[2] == 'R' && RT_flag) {
       for (i=0;i<nc_RT;i++) {
         XGRead(&ftemp, 4, 1, DMPFile, "float");
         sp_ex[i]=ftemp;
         XGRead(&ftemp, 4, 1, DMPFile, "float");
         sp_exm[i]=ftemp;
       }
       for (i=0;i<2;i++) {
         for (j=0;j<nc_RT;j++) {
           XGRead(&ftemp, 4, 1, DMPFile, "float");
           ex_rate_show[i][j]=ftemp;
         }
       }
       XGRead(&ftemp, 4, 1, DMPFile, "float");
       it_rate=(int) ftemp;
     }
   /*-------------------------------------------------*/
   fclose(DMPFile);
   printf("End Restore\n");
 }

#else
  if ((DMPFile = fopen(filename, "r+b")) == NULL)
  {
    puts("Restore: open failed");
    return;
  }
  
  XGRead(Rev, 1, 4, DMPFile, "char");
  for (i=0; i<4; i++)
    if (Rev[i]!=Revision[i] && Rev[i]!=RevisionR[i]) {
      puts("Incompatible dump file version");
      putchar(7);
      exit(1);
    }

  XGRead(&ftemp, 4, 1, DMPFile, "float");
  atime = ftemp;
  XGRead(&sigma, 4, 1, DMPFile, "float");
  XGRead(&q_0, 4, 1, DMPFile, "float");
  XGRead(&q_1, 4, 1, DMPFile, "float");
  XGRead(&q_2, 4, 1, DMPFile, "float");
  XGRead(&q_3, 4, 1, DMPFile, "float");
  XGRead(&nsp, 4, 1, DMPFile, "int");
  XGRead(iwall, 4, nsp, DMPFile, "float");
  XGRead(np, 4, nsp, DMPFile, "int");
  XGRead(&nc2p, 4, 1, DMPFile, "float");
  
  printf("Reading Restore:\n");
  printf("time = %g\n", atime);
  printf("nsp = %d\n", nsp);
  printf("np[0] = %d\n", np[0]);
  printf("np[1] = %d\n", np[1]);
     
  oldsigma=sigma;   /* This is a cludge until I can put oldsigma into dump */
  length = r1 - r0;
  for (j=0; j<nsp; j++)
  {
    for (i=0; i<np[j]; i++)
    {
      XGRead(&ftemp, 4, 1, DMPFile, "float");
      r[j][i] = length*ftemp +r0;
      XGRead(&ftemp, 4, 1, DMPFile, "float");
      v_r[j][i] = ftemp/vscale;
      XGRead(&ftemp, 4, 1, DMPFile, "float");
      v_theta[j][i] = ftemp/vscale;
      XGRead(&ftemp, 4, 1, DMPFile, "float");
      v_z[j][i] = ftemp/vscale;
    }
  }

  /*------ Load data for Radiation Transport --------*/
  if (Rev[2] == 'R' && RT_flag) {
    for (i=0;i<nc_RT;i++) {
      XGRead(&ftemp, 4, 1, DMPFile, "float");
      sp_ex[i]=ftemp;
      XGRead(&ftemp, 4, 1, DMPFile, "float");
      sp_exm[i]=ftemp;
    }
    for (i=0;i<2;i++) {
      for (j=0;j<nc_RT;j++) {
        XGRead(&ftemp, 4, 1, DMPFile, "float");
        ex_rate_show[i][j]=ftemp;
      }
    }
    XGRead(&ftemp, 4, 1, DMPFile, "float");
    it_rate=(int) ftemp;
  }

  fclose(DMPFile);
  printf("End Restore\n");
#endif

}

void Quit(void)
{
#ifdef MPI_1D
  endwtime = MPI_Wtime();
  if (my_rank == ROOT){
    printf("wall clock run time = %f\n", endwtime-time0);
    printf("init time = %f\n", time0-startwtime);
  }
  
  if (my_rank == ROOT) printf("Finalizing MPI\n");
  MPI_Finalize();
  exit(0);
#else
  endwtime = clock();
  printf("wall clock run time = %f\n", (float)(endwtime-time0)/CLOCKS_PER_SEC);
  printf("init time = %f\n", (float)(time0-startwtime)/CLOCKS_PER_SEC);
#endif
}