Add variable Np support to the fix emit/surf command

doc/fix_emit_surf.html

@@ -28,6 +28,7 @@ <H3>fix emit/surf command
 <LI>keyword = <I>n</I> or <I>normal</I> or <I>nevery</I> or <I>perspecies</I> or <I>region</I> or <I>subsonic</I> 
 
 <PRE>  <I>n</I> value = Np = number of particles to create
+                   Np can be a variable (see below)
   <I>normal</I> value = yes or no = emit normal to surface elements or with streaming velocity
   <I>nevery</I> value = Nstep = add particles every this many timesteps
   <I>perspecies</I> value = <I>yes</I> or <I>no</I>
@@ -71,9 +72,11 @@ <H3>fix emit/surf command
 within a grid cell is based on this flux and additional global, flow,
 and surface element properties:
 </P>
-<UL><LI>global property: <I>fnum</I> ratio as specified by the <A HREF = "global.html"">global</A> command
-<LI>flow properties: number density, streaming velocity, and thermal temperature
-<LI>surface element properties: portion of surface element area that overlaps with the grid cell and its orientation relative to the streaming velocity 
+<UL><LI>global property: <I>fnum</I> ratio as specified by the
+<LI><A HREF = "global.html"">global</A> command flow properties: number density,
+<LI>streaming velocity, and thermal temperature surface element
+<LI>properties: portion of surface element area that overlaps with the
+<LI>grid cell and its orientation relative to the streaming velocity 
 </UL>
 <P>The flow properties are defined for the specified mixture via the
 <A HREF = "mixture.html">mixture</A> command.
@@ -116,11 +119,28 @@ <H3>fix emit/surf command
 <P>The <I>n</I> keyword can alter how many particles are added, which can be
 useful for debugging purposes.  If <I>Np</I> is set to 0, then the number
 of added particles is a function of <I>fnum</I>, <I>nrho</I>, and other mixture
-settings, as described above.  If <I>Np</I> is set to a value > 0, then the
-<I>fnum</I> and <I>nrho</I> settings are ignored, and exactly <I>Np</I> particles are
-added on each insertion timestep.  This is done by dividing <I>Np</I> by
-the total number of grid cell/surface element pairs and adding an
-equal number of particles per pair.
+settings, as described above.
+</P>
+<P>If <I>Np</I> is set to a value > 0, then the <I>fnum</I> and <I>nrho</I> settings are
+ignored, and roughly <I>Np</I> particles are added on each insertion
+timestep.  For each grid cell/surface element pair, its target number
+of emitted particles is set to its fraction of the total emission area
+(for all grid cell/surface element pairs), multiplied by <I>Np</I>.  If
+that results in a fractional value, then an extra particle is emitted
+depending on the value of a random number, as explained above.
+</P>
+<P>The <I>Np</I> value can be also be specified as an equal-style
+<A HREF = "variable.html">variable</A>.  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case,
+the variable will be evaluated on each emission timestep, and its
+value used as <I>Np</I> on that step to determine the target number of
+emitted particles for each grid cell/surface element pair, the same as
+described in the preceeding paragraph.
+</P>
+<P>Equal-style variables can specify formulas with various mathematical
+functions, and include <A HREF = "status_style.html">stats_style</A> command
+keywords for the simulation box parameters and timestep and elapsed
+time.  Thus it is easy to specify a time-dependent value of <I>Np</I>.
 </P>
 <P>The <I>normal</I> keyword can be used to alter how velocities are set for
 added particles.  If <I>normal</I> is set to <I>no</I>, then a particle's
@@ -240,8 +260,8 @@ <H3>fix emit/surf command
 </P>
 <P><B>Restrictions:</B>
 </P>
-<P>A <I>n</I> setting of <I>Np</I> > 0 can only be used with a <I>perspecies</I> setting
-of <I>no</I>.
+<P>A <I>n</I> setting of <I>Np</I> > 0 or <I>Np</I> as a variable can only be used with
+a <I>perspecies</I> setting of <I>no</I>.
 </P>
 <P>If <I>normal</I> is set to <I>no</I>, which is the default, then unlike the <A HREF = "fix_emit/face.html">fix
 emit/face</A> command, no warning is issued if a

doc/fix_emit_surf.txt

@@ -19,6 +19,7 @@ group-ID = ID of surface group that emits particles :l
 zero or more keyword/value pairs may be appended :l
 keyword = {n} or {normal} or {nevery} or {perspecies} or {region} or {subsonic} :l
   {n} value = Np = number of particles to create
+                   Np can be a variable (see below)
   {normal} value = yes or no = emit normal to surface elements or with streaming velocity
   {nevery} value = Nstep = add particles every this many timesteps
   {perspecies} value = {yes} or {no}
@@ -62,9 +63,11 @@ given by equation 4.22 of "(Bird94)"_#Bird94.  The number of particles
 within a grid cell is based on this flux and additional global, flow,
 and surface element properties:
 
-global property: {fnum} ratio as specified by the "global"_global.html" command
-flow properties: number density, streaming velocity, and thermal temperature
-surface element properties: portion of surface element area that overlaps with the grid cell and its orientation relative to the streaming velocity :ul
+global property: {fnum} ratio as specified by the
+"global"_global.html" command flow properties: number density,
+streaming velocity, and thermal temperature surface element
+properties: portion of surface element area that overlaps with the
+grid cell and its orientation relative to the streaming velocity :ul
 
 The flow properties are defined for the specified mixture via the
 "mixture"_mixture.html command.
@@ -107,11 +110,28 @@ of particles emitting from each surface element.
 The {n} keyword can alter how many particles are added, which can be
 useful for debugging purposes.  If {Np} is set to 0, then the number
 of added particles is a function of {fnum}, {nrho}, and other mixture
-settings, as described above.  If {Np} is set to a value > 0, then the
-{fnum} and {nrho} settings are ignored, and exactly {Np} particles are
-added on each insertion timestep.  This is done by dividing {Np} by
-the total number of grid cell/surface element pairs and adding an
-equal number of particles per pair.
+settings, as described above.
+
+If {Np} is set to a value > 0, then the {fnum} and {nrho} settings are
+ignored, and roughly {Np} particles are added on each insertion
+timestep.  For each grid cell/surface element pair, its target number
+of emitted particles is set to its fraction of the total emission area
+(for all grid cell/surface element pairs), multiplied by {Np}.  If
+that results in a fractional value, then an extra particle is emitted
+depending on the value of a random number, as explained above.
+
+The {Np} value can be also be specified as an equal-style
+"variable"_variable.html.  If the value is a variable, it should be
+specified as v_name, where name is the variable name.  In this case,
+the variable will be evaluated on each emission timestep, and its
+value used as {Np} on that step to determine the target number of
+emitted particles for each grid cell/surface element pair, the same as
+described in the preceeding paragraph.
+
+Equal-style variables can specify formulas with various mathematical
+functions, and include "stats_style"_status_style.html command
+keywords for the simulation box parameters and timestep and elapsed
+time.  Thus it is easy to specify a time-dependent value of {Np}.
 
 The {normal} keyword can be used to alter how velocities are set for
 added particles.  If {normal} is set to {no}, then a particle's
@@ -231,8 +251,8 @@ a run is performed.
 
 [Restrictions:]
 
-A {n} setting of {Np} > 0 can only be used with a {perspecies} setting
-of {no}.
+A {n} setting of {Np} > 0 or {Np} as a variable can only be used with
+a {perspecies} setting of {no}.
 
 If {normal} is set to {no}, which is the default, then unlike the "fix
 emit/face"_fix_emit/face.html command, no warning is issued if a

src/create_particles.cpp

@@ -319,7 +319,7 @@ void CreateParticles::command(int narg, char **arg)
   double time1 = MPI_Wtime();
 
   bigint nprevious = particle->nglobal;
-  
+
   if (single) create_single();
   else if (!globalflag) {
     if (twopass) create_local_twopass();
@@ -459,22 +459,22 @@ void CreateParticles::create_local()
   // insertvol = subset of flowvol for cells eligible for insertion
   //   insertvol = flowvol if cutflag = 1
   //   insertvol < flowvol possible if cutflag = 0 (no cut cells)
-  
+
   double flowvolme = 0.0;
   double insertvolme = 0.0;
-  
+
   for (int icell = 0; icell < nglocal; icell++) {
     if (cinfo[icell].type == INSIDE) continue;
     if (cells[icell].nsplit > 1) continue;
     if (region && region->bboxflag &&
         outside_region(dimension,cells[icell].lo,cells[icell].hi))
       continue;
-    
+
     flowvolme += cinfo[icell].volume / cinfo[icell].weight;
     if (!cutflag && cells[icell].nsurf) continue;
     insertvolme += cinfo[icell].volume / cinfo[icell].weight;
   }
-  
+
   // calculate total Np if not set explicitly
   // based on total flowvol and mixture density
 
@@ -483,9 +483,9 @@ void CreateParticles::create_local()
     MPI_Allreduce(&flowvolme,&flowvol,1,MPI_DOUBLE,MPI_SUM,world);
     np = particle->mixture[imix]->nrho * flowvol / update->fnum;
   }
-  
+
   // gather cummulative insertion volumes across all procs
-  
+
   double volupto;
   MPI_Scan(&insertvolme,&volupto,1,MPI_DOUBLE,MPI_SUM,world);
 
@@ -542,7 +542,7 @@ void CreateParticles::create_local()
   double x[3],v[3],xcell[3],vstream_variable[3];
   double ntarget,scale,rn,vn,vr,theta1,theta2,erot,evib;
   double *lo,*hi;
-  
+
   double tempscale = 1.0;
   double sqrttempscale = 1.0;
 
@@ -557,9 +557,9 @@ void CreateParticles::create_local()
         outside_region(dimension,cells[icell].lo,cells[icell].hi))
       continue;
     if (!cutflag && cells[icell].nsurf) continue;
-    
+
     volsum += cinfo[icell].volume / cinfo[icell].weight;
-    
+
     ntarget = nme * volsum/insertvolme - nprev;
     npercell = static_cast<int> (ntarget);
 
@@ -568,7 +568,7 @@ void CreateParticles::create_local()
 
     lo = cells[icell].lo;
     hi = cells[icell].hi;
-    
+
     if (densflag) {
       scale = density_variable(lo,hi);
       ntarget *= scale;
@@ -580,7 +580,7 @@ void CreateParticles::create_local()
 
     if (cells[icell].nsurf)
       pflag = grid->point_outside_surfs(icell,xcell);
-    
+
     for (int m = 0; m < ncreate; m++) {
 
       // generate random position X for new particle
@@ -593,7 +593,7 @@ void CreateParticles::create_local()
       // if surfs, check if random position is in flow region
       // if subcell, also check if in correct subcell
       // if not, attempt new insertion up to MAXATTEMPT times
-      
+
       if (cells[icell].nsurf && pflag) {
         int nattempt = 0;
         while (nattempt < MAXATTEMPT) {
@@ -606,15 +606,15 @@ void CreateParticles::create_local()
           }
 
           nattempt++;
-          
+
           x[0] = lo[0] + random->uniform() * (hi[0]-lo[0]);
           x[1] = lo[1] + random->uniform() * (hi[1]-lo[1]);
           x[2] = lo[2] + random->uniform() * (hi[2]-lo[2]);
           if (dimension == 2) x[2] = 0.0;
         }
 
         // particle insertion was unsuccessful
-        
+
         if (nattempt >= MAXATTEMPT) continue;
       }
 
@@ -662,7 +662,7 @@ void CreateParticles::create_local()
       evib = particle->evib(ispecies,temp_vib*tempscale,random);
 
       id = MAXSMALLINT*random->uniform();
-      
+
       particle->add_particle(id,ispecies,icell,x,v,erot,evib);
 
       if (nfix_update_custom)
@@ -708,17 +708,17 @@ void CreateParticles::create_local_twopass()
   // insertvol = subset of flowvol for cells eligible for insertion
   //   insertvol = flowvol if cutflag = 1
   //   insertvol < flowvol possible if cutflag = 0 (no cut cells)
-  
+
   double flowvolme = 0.0;
   double insertvolme = 0.0;
-  
+
   for (int icell = 0; icell < nglocal; icell++) {
     if (cinfo[icell].type == INSIDE) continue;
     if (cells[icell].nsplit > 1) continue;
     if (region && region->bboxflag &&
         outside_region(dimension,cells[icell].lo,cells[icell].hi))
       continue;
-    
+
     flowvolme += cinfo[icell].volume / cinfo[icell].weight;
     if (!cutflag && cells[icell].nsurf) continue;
     insertvolme += cinfo[icell].volume / cinfo[icell].weight;
@@ -732,9 +732,9 @@ void CreateParticles::create_local_twopass()
     MPI_Allreduce(&flowvolme,&flowvol,1,MPI_DOUBLE,MPI_SUM,world);
     np = particle->mixture[imix]->nrho * flowvol / update->fnum;
   }
-  
+
   // gather cummulative insertion volumes across all procs
-  
+
   double volupto;
   MPI_Scan(&insertvolme,&volupto,1,MPI_DOUBLE,MPI_SUM,world);
 
@@ -799,7 +799,7 @@ void CreateParticles::create_local_twopass()
 
   // first pass, just calculate # of particles to create
   // ncreate_values[icell] = # of particles to create in ICELL
-  
+
   int *ncreate_values;
   memory->create(ncreate_values, nglocal, "create_particles:ncreate");
 
@@ -813,7 +813,7 @@ void CreateParticles::create_local_twopass()
     if (!cutflag && cells[icell].nsurf) continue;
 
     volsum += cinfo[icell].volume / cinfo[icell].weight;
-    
+
     ntarget = nme * volsum/insertvolme - nprev;
     npercell = static_cast<int> (ntarget);
 
@@ -839,7 +839,7 @@ void CreateParticles::create_local_twopass()
   }
 
   // second pass, create particles using ncreate_values
-  
+
   for (int icell = 0; icell < nglocal; icell++) {
     if (cinfo[icell].type == INSIDE) continue;
     if (cells[icell].nsplit > 1) continue;
@@ -858,7 +858,7 @@ void CreateParticles::create_local_twopass()
 
     if (cells[icell].nsurf)
       pflag = grid->point_outside_surfs(icell,xcell);
-    
+
     for (int m = 0; m < ncreate; m++) {
 
       // generate random position X for new particle
@@ -871,7 +871,7 @@ void CreateParticles::create_local_twopass()
       // if surfs, check if random position is in flow region
       // if subcell, also check if in correct subcell
       // if not, attempt new insertion up to MAXATTEMPT times
-      
+
       if (cells[icell].nsurf && pflag) {
         int nattempt = 0;
         while (nattempt < MAXATTEMPT) {
@@ -884,15 +884,15 @@ void CreateParticles::create_local_twopass()
           }
 
           nattempt++;
-          
+
           x[0] = lo[0] + random->uniform() * (hi[0]-lo[0]);
           x[1] = lo[1] + random->uniform() * (hi[1]-lo[1]);
           x[2] = lo[2] + random->uniform() * (hi[2]-lo[2]);
           if (dimension == 2) x[2] = 0.0;
         }
 
         // particle insertion was unsuccessful
-        
+
         if (nattempt >= MAXATTEMPT) continue;
       }
 
@@ -940,7 +940,7 @@ void CreateParticles::create_local_twopass()
       evib = particle->evib(ispecies,temp_vib*tempscale,random);
 
       id = MAXSMALLINT*random->uniform();
-      
+
       particle->add_particle(id,ispecies,icell,x,v,erot,evib);
 
       if (nfix_update_custom)
@@ -950,7 +950,7 @@ void CreateParticles::create_local_twopass()
   }
 
   // clean up
-  
+
   memory->destroy(ncreate_values);
 
   delete random;

src/cut2d.h

@@ -81,11 +81,11 @@ class Cut2d : protected Pointers {
   int split(cellint, double *, double *, int, surfint *,
             double *&, int *, int *, int &, double *);
   int split_face(int, int, double *, double *);
-  
+
   int clip_external(double *, double *, double *, double *, double *);
   int sameedge(double *, double *);
   int sameedge_external(double *, double *, double *, double *);
-  
+
  private:
   int axisymmetric;
   int implicit;
@@ -103,7 +103,7 @@ class Cut2d : protected Pointers {
   MyVec<int> used;       // 0/1 flag for each point when walking loops
 
   // methods
-  
+
   void build_clines();
   int weiler_build();
   void weiler_loops();