Merge pull request #337 from paulromano/evaporation

Improved evaporation spectrum sampling algorithm

docs/source/methods/physics.rst

@@ -682,17 +682,20 @@ nuclear temperature, which is a function of the incoming energy of the
 neutron. The ACE format contains a list of nuclear temperatures versus incoming
 energies. The nuclear temperature is interpolated between neighboring incoming
 energies using a specified interpolation law. Once the temperature :math:`T` is
-determined, we then calculate a candidate outgoing energy based on rule C45 in
-the `Monte Carlo Sampler`_:
+determined, we then calculate a candidate outgoing energy based on the algorithm
+given in LA-UR-14-27694_:
 
 .. math::
     :label: evaporation-E
 
-    E' = -T \log (\xi_1 \xi_2)
+    E' = -T \log ((1 - g\xi_1)(1 - g\xi_2))
 
-where :math:`\xi_1, \xi_2` are random numbers sampled on the unit
-interval. The outgoing energy is only accepted according to a specified
-restriction energy as in equation :eq:`maxwell-restriction`.
+where :math:`g = 1 - e^{-w}`, :math:`w = (E - U)/T`, :math:`U` is the
+restriction energy, and :math:`\xi_1, \xi_2` are random numbers sampled on the
+unit interval. The outgoing energy is only accepted according to the restriction
+energy as in equation :eq:`maxwell-restriction`. This algorithm has a much
+higher rejection efficiency than the standard technique, i.e. rule C45 in the
+`Monte Carlo Sampler`_.
 
 ACE Law 11 - Energy-Dependent Watt Spectrum
 +++++++++++++++++++++++++++++++++++++++++++
@@ -1591,6 +1594,8 @@ References
 
 .. _Monte Carlo Sampler: https://laws.lanl.gov/vhosts/mcnp.lanl.gov/pdf_files/la-9721_3rdmcsampler.pdf
 
+.. _LA-UR-14-27694: http://permalink.lanl.gov/object/tr?what=info:lanl-repo/lareport/LA-UR-14-27694
+
 .. _MC21: http://www.osti.gov/bridge/servlets/purl/903083-HT5p1o/903083.pdf
 
 .. _Sutton and Brown: http://www.osti.gov/bridge/product.biblio.jsp?osti_id=307911

src/physics.F90

@@ -1832,14 +1832,15 @@ contains
       lc = 2 + 2*NR + 2*NE
       U = edist % data(lc + 1)
 
+      y = (E_in - U)/T
+      v = 1 - exp(-y)
+
       ! sample outgoing energy based on evaporation spectrum probability
       ! density function
       n_sample = 0
       do
-        r1 = prn()
-        r2 = prn()
-        E_out = -T * log(r1*r2)
-        if (E_out <= E_in - U) exit
+        x = -log((1 - v*prn())*(1 - v*prn()))
+        if (x <= y) exit
 
         ! check for large number of rejections
         n_sample = n_sample + 1
@@ -1849,6 +1850,8 @@ contains
         end if
       end do
 
+      E_out = x*T
+
     case (11)
       ! =======================================================================
       ! ENERGY-DEPENDENT WATT SPECTRUM