start.c

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

  }

  r_tot  = (float **)malloc(nsp*sizeof(float *));
  vr_tot = (float **)malloc(nsp*sizeof(float *));
  vtheta_tot = (float **)malloc(nsp*sizeof(float *));
  vz_tot = (float **)malloc(nsp*sizeof(float *));

  for(isp=0; isp<nsp; isp++) {
    r_tot[isp]  = (float *)malloc(num_procs*maxnp[isp]*sizeof(float));
    vr_tot[isp] = (float *)malloc(num_procs*maxnp[isp]*sizeof(float));
    vtheta_tot[isp] = (float *)malloc(num_procs*maxnp[isp]*sizeof(float));
    vz_tot[isp] = (float *)malloc(num_procs*maxnp[isp]*sizeof(float));

    for(j=0; j<maxnp[isp]; j++)
      r_tot[isp][j]=vr_tot[isp][j]=vtheta_tot[isp][j]=vz_tot[isp][j]=0.0;

  }  

  if (nsp && (!r_tot[nsp-1] || !vr_tot[nsp-1] || !vtheta_tot[nsp-1] ||
              !vz_tot[nsp-1])) {
    if(my_rank == ROOT)
      puts("START: r, v malloc failed");
    exit(-1);
  } 

#endif /* MPI_1D */
  

  /*********************************************************/
  /* Load in all species' particles, properly distributed. */
  /* Start from a dump file if given otherwise load the    */
  /* system with the initial distribution                  */

  if (WasDumpFileGiven) 
  {
       ntemp = nc2p;
       Restore(theDumpFile);

 /* Rescale normalisation variables using nc2p from the Dump file */
       for (isp=0; isp<nsp; isp++) intoi[isp] *= nc2p/ntemp; 
       rhon *= nc2p/ntemp; 
  }
  else if (profile == 1)
       unif_load(initnp);
  else if (profile == 2)
     bes_load(initnp,iinitn);
  else if (profile == 3)
     cos_load(initnp,iinitn);


}	/* End START */

/***************************************************************/
/*   If species parameters exist read them from input file    
**   Shape of density profile determined by parameter profile:
**   profile = 1   uniform profile
**   profile = 2   bessel function profile
**   profile = 3   cosine profile
*/

void species(int initnp[NSMAX][2], int isp)
{
  char aachar[200];
  int np_start, ii;
  float j0_l, j0_r, v0r_l, v0r_r, vcr_l, vcr_r, vtr_l, vtr_r,
        v0_t, vt_t, v0z, vtz;
  float c0, temin, temax, temin_mid, temax_mid, area1;
  
  /* Read in SPECIES parameters... */
  while (fscanf(InputDeck, "%f %f %f %f %f %d %d", &charge[isp], &mass[isp],
		&j0_l, &j0_r, &iinitn, &profile, &sp_k[isp]) <7)
    fscanf(InputDeck, "%s", aachar);
  
  while (fscanf(InputDeck, "%f %f %f %f %f %f ", &v0r_l, &v0r_r, &vtr_l,
		&vtr_r, &vcr_l, &vcr_r) <6)
    fscanf(InputDeck, "%s", aachar);
  
  while (fscanf(InputDeck, "%f %f %f %f", &v0_t, &vt_t, &v0z, &vtz) <4)
    fscanf(InputDeck, "%s", aachar);

/* Energy distribution at outside wall */
  while (fscanf(InputDeck, "%d %f %f %d", &nbin[isp], &temin,
		&temax, &maxnp[isp]) <4)
    fscanf(InputDeck, "%s", aachar);

/* Energy distribution at plasma centre */
  while (fscanf(InputDeck, "%d %g %g %g %g", &nbin_mid[isp], &temin_mid,
                                           &temax_mid, &rs_mid[isp], &rf_mid[isp]) <5)
     fscanf(InputDeck, "%s", aachar);

/* Velocity distribution at plasma centre */
  while (fscanf(InputDeck, "%g %g %d %g %g %d %g %g %d", &vr_l[isp], &vr_u[isp], &nvrbin[isp],
            &vt_l[isp], &vt_u[isp], &nvtbin[isp], &vz_l[isp], &vz_u[isp], &nvzbin[isp]) <9)
     fscanf(InputDeck, "%s", aachar);

  if (nvrbin[isp]||nvtbin[isp]||nvzbin[isp]) vel_dist_accum=1;
  
  /* Assign input parameters to arguments */
  v0[isp][0] = fabs(v0r_l);          /* drift velocity v > 0 (moving right) */
  v0[isp][1] = -fabs(v0r_r);         /* drift velocity v < 0 (moving left) */
  vc[isp][0] = fabs(vcr_l);          /* cutoff velocity for v > 0, thermal dist */
  vc[isp][1] = -fabs(vcr_r);         /* cutoff velocity for v < 0, thermal dist */
  vt[isp][0] = fabs(vtr_l);          /* thermal velocity, v > 0 */
  vt[isp][1] = -fabs(vtr_r);         /* thermal velocity, v < 0 */
  v0_th[isp]  = v0_t;                /* drift velocity in theta direction */
  vt_th[isp]  = fabs(vt_t);          /* thermal velocity in theta direction */
  v0_z[isp]  = v0z;                  /* drift velocity in z direction */
  vt_z[isp]  = fabs(vtz);            /* thermal velocity in z direction */

/* parameters for energy and velocity distribution diagnostics */
  emin[isp]  = temin;              /* min energy measured in energy dist diag */ 
  de[isp]	  = (temax-temin)/(nbin[isp]-1);      /* energy bin size */
  dtheta[isp]= 0.5*PI/(nbin[isp]-1);                    /* angular bin size */
  emin_mid[isp] = temin_mid;
  de_mid[isp] = (temax_mid-temin_mid)/(nbin_mid[isp]);
  
  /*determine max vel of right & left moving particles for making distribution function */
  vrscale= max(vrscale, v0[isp][0]+vt[isp][0]);
  vrscale= max(vrscale, -v0[isp][1]+vt[isp][1]);

  /* Check that nc2p and iinitn produce # sim particles w/in range */
  if (!WasDumpFileGiven) 
  {
     if (profile == 1) 
     {
       np_start = iinitn*reactor_vol/nc2p + 0.5;
      #ifdef MPI_1D
         if (my_rank == ROOT) 
         {   
           printf(" den_init = %e vol = %f \n",iinitn, reactor_vol);
           printf("np = %d maxnp = %d ",np_start,maxnp[isp]);
           printf("nc2p = %e \n",nc2p);
         }
      #else
         printf(" den_init = %e vol = %f \n",iinitn, reactor_vol);
         printf("np = %d maxnp = %d ",np_start,maxnp[isp]);
         printf("nc2p = %e \n",nc2p);
      #endif
     }
     else 
     {
       np_start = 0;
       c0 = iinitn/dntod;
       if (profile == 2)
       {
         for (ii=0; ii<=nc-2; ii++)
           np_start += c0*J_0((r_grid[ii]-r0)/(r1-r0))*(r_grid[ii+1]*r_grid[ii+1] - r_grid[ii]*r_grid[ii]);
       }
       else if (profile == 3)
       {
         for (ii=0; ii<=nc-2; ii++)
           np_start += c0*cos(PI*0.5*(float)(r_grid[ii]-r0)/(r1-r0))*(r_grid[ii+1]*r_grid[ii+1] - r_grid[ii]*r_grid[ii]);
       }
      #ifdef MPI_1D
        if (my_rank == ROOT) 
        {
          printf(" n_pk = %e \n",iinitn);
          printf("np = %d maxnp = %d ",np_start,maxnp[isp]);
          printf("nc2p = %e \n",nc2p);
	}
      #else
        printf(" n_pk = %e \n",iinitn);
        printf("np = %d maxnp = %d ",np_start,maxnp[isp]);
        printf("nc2p = %e \n",nc2p);
      #endif
     }   /* end of if (profile) loop */

  /* increase nc2p & reduce # simulation particles */
     if (np_start > maxnp[isp])                  
     {
         nc2p *= np_start/(0.9*maxnp[isp]);
         np_start = 0.9*maxnp[isp];
       #ifdef MPI_1D
         if (my_rank == ROOT)
           printf("SPECIES - too many simulation particles, decreasing to %d and increasing nc2p to %e  \n\n",np_start, nc2p);
       #else
         printf("SPECIES - too many simulation particles, decreasing to %d and increasing nc2p to %e  \n\n",np_start, nc2p);  
       #endif
     }
  /* decrease nc2p & increase # simulation particles*/
     if (iinitn && np_start < min_np)
     { 
         nc2p *= np_start/(1.2*min_np);
         np_start = 1.2*min_np;
      #ifdef MPI_1D
        if (my_rank == ROOT) 
          printf("SPECIES - too few simulation particles, increasing to %d and decreasing nc2p to %e \n\n",np_start, nc2p);
      #else
        printf("SPECIES - too few simulation particles, increasing to %d and decreasing nc2p to %e \n\n",np_start, nc2p);
      #endif
     }
 /* If particles are moving both left and right */
     if((vt[isp][0] || v0[isp][0]) && (vt[isp][1] || v0[isp][1])) 
     {
       initnp[isp][0] = 0.5*np_start;
       initnp[isp][1] = 0.5*np_start;
     }
/* If particles are moving only to the left */
     else if((vt[isp][1] || v0[isp][1]) && (!vt[isp][0]  && !v0[isp][0])) 
     {
       initnp[isp][0] = 0.0;
       initnp[isp][1] = np_start;
     }
/* If particles are moving only to the right, or are stationary */
     else 
     {
       initnp[isp][0] = np_start;
       initnp[isp][1] = 0.0;
     }

  }    /* end of loop for initialising density */

/* MPI: Divide particles among procs */
     #ifdef MPI_1D
       maxnp[isp] = floor(0.5 + ((float) maxnp[isp])/num_procs);
       iinitn = floor(0.5 + ((float)iinitn)/num_procs);
       initnp[isp][0] = floor(0.5 + ((float)initnp[isp][0])/num_procs);
       initnp[isp][1] = floor(0.5 + ((float)initnp[isp][1])/num_procs);
     #endif

  /* j0_l & j0_r = injected current density from left & right electrodes */
  /* Set enter = no. of simulation particles injected PER DT */
  area1 = TWOPI*r1*height;               /* area outer electrode */
  enter[isp][0] = 0.0;
  if (r0>0.0) enter[isp][0] = fabs(j0_l*area0*dt*sp_k[isp]/(charge[isp]*nc2p));
  enter[isp][1] = fabs(j0_r*area1*dt*sp_k[isp]/(charge[isp]*nc2p));

  #ifdef MPI_1D
    /* MPI: Divide injected particles among procs */
    enter[isp][0] /= num_procs;
    enter[isp][1] /= num_procs;
  #endif

  if (enter[isp][0]> 0.) inject[isp] = 1;
  if (enter[isp][1]> 0.) inject[isp] = 1;

}	/* End SPECIES */

/***************************************************************/
/* Bessel's function calc, for loading particles               */
float J_0(float z)
{
   float x, j0;

   z *= 2.4;
   x = 0.25*z*z;
   j0 = 1 - x + 0.25*x*x - x*x*x/36.0;
   return (j0);
}
/*------------------------- LEEHJ -----------------------------*/
void input_for_rad()
{
  char aachar[80];
  float k0_input;

  while (fscanf(InputDeck, "%d %g %c %g %g %d", &nc_RT, &lambda0, &Lineshape,
         &A_ki, &k0_input, &sp_k_ex)<6)
    fscanf(InputDeck, "%s", aachar);

  while (fscanf(InputDeck, "%g %g %d", &xmin, &xmax, &xbin) <3)
    fscanf(InputDeck, "%s", aachar);

  if      (Lineshape=='L' || Lineshape=='l') { lineshape=1; Cx = 1/M_PI;       }
  else if (Lineshape=='D' || Lineshape=='d') { lineshape=2; Cx = 1/sqrt(M_PI); }
  else if (Lineshape=='V' || Lineshape=='v') { lineshape=3; Cx = 1/M_PI;       }
  else {
    puts("Select proper lineshpae\n");
    exit(1);
  }
  tau=1e-8/A_ki;
  A_ki*=1e8;

  if (lineshape==1) {
    delta_nu=1.0555e-7*1.532*0.2214*lambda0*1e-9*gas_den;
    k0=lambda0*lambda0*A_ki*Cx*3/(4*M_PI*delta_nu);
    k0*=1e-18*gas_den;       /* 1e-18 = 1e-9[m]*1e-9[m]         */
  }
  else if (lineshape==2) {
    delta_nu=2*sqrt(2*log(2)*T_gas*1.602/(mass[1]*1e19))/lambda0*1e9;
    k0=lambda0*lambda0*A_ki*Cx*3/(4*M_PI*delta_nu)*sqrt(log(2));
    k0*=1e-18*gas_den;       /* 1e-18 = 1e-9[m]*1e-9[m]         */
  }
  I_factor=2/pow(lambda0,3)*6.6261*2.99792458*10/gas_den/3.0;
  lambda0*=1e-9;        /* lambda0[nm] is normalized as nm */
  printf("gas_den=%g,\t",gas_den);
  printf("delta_nu=%g GHz\t k0=%e\n",delta_nu*1e-9,k0);
  k0=k0_input;
}
/***************************************************************/
/* Routines to calculate nearest mesh point to particle position
** used in field and move. Depends on the type of mesh used in the simulation
*/
/* ------------------------------------------------------------------- */
int unif_grid(int isp, int ii)
{
   int j;
   j = r[isp][ii] - r0;
   return(j);
}
/* ------------------------------------------------------------------- */
int const_vol_grid(int isp, int ii)
{
   int j;
   j = (r[isp][ii]*r[isp][ii] - r0*r0)/(2.0*r0 + 1.0);
   return(j);
}
/* ------------------------------------------------------------------- */
int linear_grid(int isp, int ii)
{
   int j;

/* Get approximate position of nearest grid to particle at r */
   j = r[isp][ii] - r0;

/* Increment j until nearest grid point found */
   while (r_grid[j+1] < r[isp][ii] && j < ng) j++;

   return(j);
}
/* ------------------------------------------------------------------- */
int split_grid(int isp, int ii)
{
   int j;
   if (r[isp][ii] <= r_grid[nc0]) j = r[isp][ii] - r0;
   else j = nc0 + (r[isp][ii] - r_grid[nc0])*dr/dr1;
   return(j);
}
/***************************************************************/